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NOMENCLATURE 

a  Channel  half-height 

A  Cross-sectional  area  of  the  tube 

D  Tube  diameter 

g  x-oomponent  of  the  acceleration  of  gravity 

p  Fluid  pressure 

p.  Pressure  at  channel  mouth 

p  Pressure  at  far  upstream  where  fluid  velocity  Is  zero 

Ap      Total  pressure  drop  in  the  entrance  region 

(P-P0) 
P       Dimensionless  pressure,  5 

H       Dimensionless  radial  coordinate,  — - 

.-  r0 

He      Reynolds  number, 

Ren     Reynolds  number, 

/urQ 
Re      r.eynolds  number,     ~JZ~ 

rQ  Tube  radius 

t  Time 

u,v,w  x,y,z-components  of  the  velocity  field 

u  Average  velocity 

u.  Center  line  velocity 

U       Dimensionless  velocity  In  x-directlon,  -*- 

u 
uf      Velocity  at  fully  developed  condition 

u*      Deviation  velocity 

x,y,z   Rectangular  coordinates 

X       Dimensionless  x-coordinate,   — ~—  =  *S; — '- 

/a2S     Hea 

x       Entrance  length 


IV 


(xe/a) 


X  Dlmensionless  entrance  length, 

X  Dlmensionless  y-coordlnate,   J 

1/  Kinematic  viscosity 

f  Fluid  density 

M.  Dynamic  viscosity 

S  Boundary  layer  thickness 
Shear  stress  at  the  wall 


Tw 

£       Dlmensionless  x-coordlnate,   (— ~— ) 

a  u 
c*      Dlmensionless  x-coordlnate  at  matching  section 


CHAPTER  1.   INTRODUCTION 

When  a  viscous  fluid  enters  a  conduit  (a  circular  tube  or  a 
flat  duct)  the  velocity  profile  of  the  flow  changes  as  the  f low - 
progresses  from  the  entrance  to  the  interior  of  the  conduit  be- 
cause the  wall  of  the  conduit  tends  to  retard  the  flow.   Even- 
tually, the  velocity  profile  will  assume  a  form  which  is  invari- 
ant with  respect  to  the  direction  of  flow.   The  flow  in  this 
region  of  changing  velocity  profile  is  the  so-called  hydrodynamlc 
entrance  region  flow. 

A  study  of  the  velocity  field  and  friction  factor  in  the 
entrance  region  of  tubes  and  ducts  Is  of  practical  importance  and 
has  been  a  subject  of  investigation  for  the  last  hundred  years. 
In  the  study  of  the  entrance  region  flow  we  usually  wish  to  ob- 
tain expressions  for  (a)  the  velocity  distribution  at  any  sec- 
tion, (b)  the  pressure  drop  between  any  two  sections,  and  (c)  the 
value  of  the  entrance  length  at  which  the  fully  developed  flow 
may  be  said  to  be  attained. 

Although  the  governing  differential  equations  for  the  en- 
trance region  flow  can  be  easily  written,  their  exact  analytical 
solutions  have  not  as  yet  been  found.   The  celebrated  boundary 
layer  theory  was  published  in  190*1  by  Prandtl.   The  equation  of 
motion  for  the  entrance  region  flow  can  be  simplified  consider- 
ably if  such  flow  Is  approximated  by  a  model  based  on  the  bound- 
ary layer  theory.   Although  the  boundary  layer  equations  are 
much  simpler  than  the  original  equations  of  motion,  we  are  still 
unable  to  obtain  their  exact  analytical  solutions.   Thus  various 


approximation  methods  of  solution  have  been  employed  to  solve  the 
entrance  region  flow  problem.   A  thorough  review  and  classifica- 
tion of  the  literature  of  this  subject  has  been  made  by  Hwang  and 
Fan  [l]*. 

The  approximate  solutions  may  be  classified  into  four  gener- 
al categories  —  the  momentum  Integral  method  of  Schiller  [U], 
the  linearization  method  of  Langhaar  15]  and  Targ  [6  J,  the  match- 
ing method  of  Schllchtlng  [7],  and  the  finite  difference  method 
{&].      These  approximate  methods  are  explained  in  detail  in  Li's 
report  \9].      However,  the  numerical  calculations  of  the  solutions 
are  not  carried  out  in  Li's  report. 

The  purpose  of  this  report  is  the  detailed  study  of  the  en- 
trance region  flow  phenomena  by  obtaining  the  numerical  results 
with  a  high  speed  computer.   Three  methods  including  Schiller's 
momentum  integral  method  [*»•],  Langhaar's  linearizing  method  [5], 
and  Schlichting's  matching  method  [7]  are  studied  numerically. 
A  series  of  numerical  analyses  Including  the  evaluation  of  modi- 
fied Bessel  functions  [  10 J,  numerical  integration  [  11,  13.  14] 
series  solution  of  differential  equations,  eigenvalue  solution 
I 12]  and  others  are  employed  in  this  report. 

The  details  of  the  numerical  analysis,  computational  proce- 
dures, flow  charts,  FORTRAN  programs  and  numerical  results  are 
presented. 


•  Numerals  in  brackets  refer  to  references  at  the  end  of  each 
chapter. 
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CHAPTER  2.   SCHILLER'S  METHOD 

2.1.  Description  of  the  Method 

Schiller's  method  Is  also  known  as  the  Karman-Pohlhausen 
momentum-Integral  method.   When  a  viscous  fluid  enters  a  olroular 
tube  with  a  uniform  velocity  profile  across  the  cross-section,  a 
thin  boundary  layer  is  often  considered  to  have  been  formed. 
This  Is  mentioned  in  the  previous  chapter.   The  thickness  of  this 
layer  increases  with  the  distance  from  the  entry  until  It  Is 
equal  to  the  radius  of  the  pipe.   Outside  the  boundary  layer,  the 
flow  is  considered  to  be  frictionless  and  accelerated.   In  this 
method  the  velocity  profile  is  assumed  to  be  a  curve  built-up  by 
two  parabolas  and  a  straight  line  (see  Fig.  2.1);  the  vertices 
of  the  parabolas  lie  on  the  border  of  the  boundary  layer.   The 
parabolic  velocity  in  the  boundary  layer  Is  assumed  to  be 

£■  =   2(f)  -  (f)2  (1) 

u0        6      6 

where 

y  =  rQ  -  r 

and  uQ  Is  the  velocity  of  the  central  core  and  is  a  function  of 
xs  6  is  the  thickness  of  the  boundary  layer. 

The  equation  of  continuity  and  the  equation  of  motion  (the 
momentum  equation)  give  two  relations  among  u0>  x,  and  6. 
Elimination  of  6  will  lead  to  a  relation  between  uQ  and  x. 

The  momentum  theory  Is  used  here  for  the  x  direction  to  ob- 
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tain  a  momentum  Integral  equation  representing  the  equation  of 
motion.   The  control  surface  Is  formed  by  two  planes  a-b  and  c-d, 
spaced  apart  from  each  other  by  the  small  distance  dx  and  by  the 
wall  (see  Fig.  2.2).   By  means  of  the  momentum  balance  on  this 
control  surface  the  following  momentum  Integral  equation  Is  ob- 
tained. 

£j7u2(r0-y)dy]  -  «„&[/'  u(rQ  -  y)dy] 

=  -2«Vt}  *  \  ro2 12  ^  "  (t)2 1  u°  ^  '     (2) 

Substituting  the  assumed  velocity  profile  expressed  by 
equation  (1)  into  the  integral  form  of  the  continuity  equation 
and  carrying  out  the  integration,  we  obtain  a  relation  between 
the  boundary  layer  thickness,  6.  and  the  center  line  velocity, 
u0,  as  shown  below j 


-  /» - 6<i  -  h 


r0  '        V  "o 


(3) 


where  u  is  the  average  velocity. 

Substituting  equations  (1)  and  (3)  into  equation  (2),  and 
carrying  out  the  integration,  we  obtain  a  relation  between  uQ 
and  x  as  follows: 
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_  16  /OS  .  22J1   sm"1  (Stl) 
5  /l+1    10  sln   l  3  ' 

-  2Z£  sin"1  (})+$]  (*) 


where 


«l-  ^-1. 

u 

Therefore  the  velocity  profile  in  the  entrance  region  is  defined 
by  equations  (1),  (3),  and  CO.   Given  a  value  of  uQ,  the  values 
of  x  and  6  are  known  from  equations  CO  and  (3).   The  velocity 
distribution  at  the  corresponding  section  is  defined  by  equation 
(1).   The  flow  is  fully  developed  when  6  =  rQ  and  uQ  =  2u. 
Hence  l\  =  1.   Then  from  equation  (4),  the  entrance  length  is 
found  to  be 

1   =  0.02877  DReD  . 

Differentiating  both  sides  of  equation  CO  with  respect  to 

d»l  du0 

x,  we  see  that  ¥*,  that  is  ■—-,    does  not  vanish  at  l\  =  1.   This 

means  that  the  center  line  velocity  does  not  asymptotically 
approach  its  final  value,  2u,  at  the  fully  developed  section. 

For  the  calculation  of  the  pressure  drop  in  the  entrance  re- 
gion, the  Bernoulli  equation  is  applied  since  it  is  assumed  that 
a  potential  flow  In  the  central  core  exists  throughout  the  en- 
trance region.   Therefore,  the  pressure  distribution  is  assumed 
to  be  dependent  only  on  the  potential  flow.   The  pressure  drop 
between  any  two  locations  in  this  region  is  expressed  as 


10 


Pl   -   P2 


fl(u0)22-  (uo)l2]  =  i&[e\z2-\z) 

+  Z{\  -\)]  •  (5) 


Since  i  =  0  andVl  =  0  at  the  entry,  and  x  ■  x  and  1^=  1  at  the 
fully  developed  section,  the  total  pressure  drop  In  the  entrance 
region  Is 

AP   =   ff  U2  . 

The  pressure  drop  In  the  entrance  region  Is  larger  than  that 
In  the  fully  developed  region  of  equal  length  because  of  the  In- 
creased friction  loss  and  the  fact  that  the  kinetic  energy  of  the 
fluid  Increases  as  it  passes  downstream.   For  the  evaluation  of 
total  pressure  drop  between  a  section  in  the  entrance  region  and 
a  section  in  the  fully  developed  region.  It  is  more  convenient  to 
treat  the  entrance  region  as  a  fully  developed  one  and  the  pres- 
sure drop  thus  calculated  is  modified  by  a  correction  factor. 
The  equation  for  calculating  the  total  pressure  drop  in  the  fully 
developed  section  Is  found  to  be 

p0  I  p      64  x   .  ., 
-££-  =  Re^D*1'16 

2 

where  1.16  is  the  correction  factor.   If  the  fluid  enters  the 
pipe  from  a  region  where  the  pressure  is  p  and  the  velocity  is 

CO 

negligible,  there  Is  a  pressure  drop  of  yf u  ,  and  then  the 
total  pressure  drop  is 


11 


^^+  2.16 


2 

The  derivation  and  solution  are  given  in  detail  in  Li's  report 
12]. 


12 

2.2.  Computational  Procedures 

To  normalize  the  equations  (1),  (3),  CO  and  (5).  the 
following  dlmenslonless  parameters  are  introduced. 


Uc 
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u 
u 


Y   _   X. 

r0 


D  .   =  J- 
el     rQ 


fuD  =   HeD 


'       /u2       ' 
Equations   (1),    (3),    CO    and   (5)    become 

f   -    ^BT*  "  (dT)2  '  (6) 

o  el  el 


D 
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x     -     2  -  y4  -  6(1  -  J-j    ,  (7) 


-  X£  sin"1   (^i)    -  3gfE  sin"1    (i)    ♦  f  ]  .  (8) 

P     -     |[^2   -Iq)    ♦   2(1  -V]   '  (9) 
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where 

*\-  B0  -  1  .  (10) 

In  order  to  use  available  subroutines,  arcslne  functions  are 
changed  to  arotangent  functions  as  follows: 

al,,"1  (2VO.,       .    -1        21-1 

sin   ( — — )  m     tan   ( -)  . 

2  J ^2+n-a2 

sin-1  (i)   =  tan-1  (-M  . 
J  2-/2 

The  computational  procedure  is 

(I)  Given  a  value  of  U„,  the  center  line  velocity, 

c 

then  \   is  obtained  from  equation  (10). 

(II)  The  values  of  the  boundary  layer  thickness,  D -, 
the  location  from  the  entry,  X,  and  the  pressure 
drop,  P,  can  be  obtained  from  equations  (7),  (8) 
and  (9)t  respectively. 

(ill)  The  velocity  distribution,  U,  at  the  corresponding 
section  Is  obtained  from  equation  (6). 

The  center  line  velocities  (or  core  velocities)  are  U  =  1.0 
at  the  entry  where  X  =  0  and  U  =  2.0  at  the  fully  developed 
section.   The  mesh  size  for  the  calculation  of  the  velocity  pro- 
file and  pressure  drop  Is  shown  in  Table  2.1. 

A  flow  chart  (Table  2.2)  will  be  followed  by  a  symbol  table 
(Table  2.3)  and  FORTRAN  program  (Table  2.^).  The  computation  Is 
carried  out  by  use  of  an  IBM  1620  computer. 


Table  2.1.   Mesh  sizes  for  the  velocity  profile  and 
pressure  drop  in  the  entrance  region 
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Uc 

*Uc 

AY 

N+ 

1.0005  , 

.0005 

Del/5 

5 

1.01 

.001 

Del/5 

5 

1.02 
1.07 
1.20 
2.0 

.005 

Del/5 

5 

.01 

Del/5 

5 

.05 

Del/10 

10 

(fully  developed) 

♦  N  =  points  within  the  boundary  layer. 
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Table  2.2.   Flow  chart  of  Schiller's    Method 
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Table   2.3.      Symbol  Table 
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FORTRAN  Symbol 

Variables 

DEL 

el 

DLTYB 

A  Y 

ET 

n. 

P 

p 

3 

& 

U 

u 

UC 

Uo 

XD 

X 

YR 

Y 
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61 
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FCRMATI/E12.5.F12.4,    itETl?.  5  ) 

FORMAT     (4RX,    2F12.5     ) 

UC«1.1  005 

5=1.4142135 

ET  =  UC-1.'- 

XD=l./16.«(58./lb.*ET-22./5.*LClJ<  l.+ET)  ) 

Xr)=X0-17./(16.*lb.)*buRT(4.  +  2.*LT-^.*I:T*tT) 

X0=XD-1./5.*5URT<  <4.-2.*£T  )/(  l.+LT)  ) 

XD=XD-1./16.*3.7*S*ATAN(  (£T-.5)/S(jKT(2.  +  ET-tT*t-T)  ) 

XD=XD-1./16.«(3.7*S*ATAN< l./(2.*S) I-2  6./3. ) 

p=(ET*ET+2.*ET)/2.. 

DEL=2.-50RT(-2.+6./UC) 

IFIUC-1. 2  I     4,4,5 

DLTYR=DEL/5. 

GC    TC    4  1 

DLTYR=DEL/li  . 

YR'DLTYR 

U=UC*(2.*(YR/DEL)-( YR/DEL) *( YR/DEL) ) 

IF(YR-OEL)  61,  61,  6 

1F( YR-ULTYR)  62,62,63 

PUNCh  1,  uC.XD.P.tjLL.YR,  o 

GC  TC  7 

PUNCH  1  ..),  YR,U 

YR«YR+DLTYR 

GC  TC  6 

IF  [UC-] 

UC=UC+.l 

GC  TC  2 

IFIUC-1. 

UC=UC+.u 

GC  TC  2 

IF  (Uf-1 

UC=UC+.0 

G-:  TC  2 

IFIUC-1 .2 

UC=UC+.01 

GC  TC  ? 

IFIUC-2.0) 

L'C  =  UC  +  .05 

GC  TC  2 

STCP 

END 
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2.3.  Results 

Some  numerical  results  are  tabulated  In  Table  2.5*   The 
developing  velocity  profile  is  shown  in  Pig.  2. J.      The  pressure 
drops  in  the  entrance  region  are  plotted  in  Pig.  2.^. 
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TABLE  ?t5.   SAMPLE  NUMERICAL  RESULTS  FOR  SCHILLER'S  METHCD 


UC  XD  P  DEL  YR  U 

.10500E+01   U.5G57E-04  0.51250E-01  0.72752E-01  Q.14550E-03  0.37800 

•  29101E-01  CJ. 67200 

.43651E-01  0.88200 

.58201E-O1  0.10080E+01 

.72752E-01  0.10500E+01 

.1100CE+01   0.2146E-03  0.10500      0.14136      0.28272E-01  0.39600 

.56544E-01  0.70400 

•84816E-01  C. 92400 

.11309  0.10560E+01 

.14136  0.11000E  +  01 

.120OOE+O1   U.9393E-03  0.22000      0.2679b      Q.&3590E-O1  0. 43200 

.10718  0.76800 

.16077  0.1008CE+01 

. 21436  0.11520E+01 

.26795  0.12000E+01 

.1300CE+01   U.2251E-02  0.34500      0.38279      G.38279E-01  0.247nn 

.76557E-01  0.468CC 

.114  8  4  0.66300 

.15311  C.832C0 

.19139  0.97500 

.22967  0.1092CE+01 

.26795  0.11830E+01 

.30623  0.12480E+01 

.34451  0.12870E+01 

.38279  0.13000E+C1 

.14C0OE  +  -1   (..4185E-C2  0.48000      0.48814      0.48814E-0]  0.2660C 

•97628E-01  0.5040C 
.14644      C. 71400 
.19526      0.89600 
.24407      0.105C0E+01 
.29289      0.11760E+01 
.34170      0.12740E+01 
.39051      0.13440E+01 
.43933      0.1386CE+01 
.48814      0.14C00E  +  01- 

•  151-0CF  +  01   W.6755E-02  0.62500      0.50579      0.58579E-01  0.28500 

.11716  0.54000 

.17574  0.76500 

.23431  0.96CC0 

.29289  G.11250E+01 

.35147  0.12600E+01 

.41005  0.13650E+01 

.46863  0.14400E+01 


.52721      0.14850E+01 
.58579      0.15000E+01 
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2.4.  Sample  Calculation  with  Dimensions 

As  soon  as  the  dimenslonless  results  are  obtained,  the 
values  with  any  dimension  can  be  calculated. 

For  example,  we  have  three  different  kinds  of  fluids,  air, 
water,  and  lubricating  oil,  with  their  properties  listed  In  Table 
2.6.  The  diameter  of  the  tube  and  the  average  velocity  are  as- 
sumed to  be  the  same  for  all  three  fluids.   Assuming  that  D  =  1/4 
ft,  u  =  1/10  ft/sec,  we  try  to  compare  the  entrance  length  and 
pressure  drop  from  the  entry  to  the  place  where  the  fluid  Is 
fully  developed.   Here,'  the  assumption  Is  made  that  the  fluids 
are  Incompressible. 

For  air: 

ReD    =    TT    -    ioTT  1.460  '  10"    =    1'712  "  1q2 

Using  the  output  data  that  X  =  .02877,  and  P  =  1.5,  we  obtain 


x       =     X   •    Hen   •   D 
e  D 


.02877   •    1.712   .   102   .   £     =     1.231  ft   , 


xe/D     =     1.231/(1/4)      =     4.924   , 
Ap     =     pQ   -  Pe     =     P  •  f  u 

=     1.5   •    (.00237)    •'  (^)2   •  jfa; 

=     2.464   .    10"7  psla 

The  calculations  for  water  and  oil  proceed  In  the  same  way, 
and  the  results  are  shown  In  Table  2.6. 
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CHAPTER  3.   LANGHAAR* S  LINEARIZING 
APPROXIMATE  METHOD 

3.1.  Description  of  the  Method  [l] 

Because  of  the  nonlinear  convectlve  terms  Involved  in  the 
boundary  layer  equations,  it  is  difficult  to  obtain  an  exact 
analytical  solution.   By  means  of  a  linearizing  approximation, 
Langhaar  assumed  that  the  convectlve  terms  in  the  boundary  layer 
equation  in  the  axial  direction,  x,  (Fig.  2.1)  are  equal  to  the 
product  of  the  kinematic  viscosity  V,  the  axial  velocity  u  and  a 
parameter  $  which  is  considered  to  be  a  function  of  x  only. 

Langhaar' s  assumption  can  be  written  as 

3u     au     3u        2  / . , 

U?I+  V^+  wii  "  **   u  (1) 

and  gives  rise  to  the  following  linearized  form  of  the  boundary 
layer  equation 

"e2u    =    sx  -7g+*^*^!>  •  (2) 

x      f  ax  ay2   3z2 

In  order  to  solve  the  problem  for  a  flow  in  a  circular  tube, 
equation  (2)  is  transformed  into  a  cylindrical  coordinate  form  as 

d2u   1  du   02  ,„, 

7" 2  +  r   d?  "  *  u  =  *  (3) 

dr 

where 


«   *Bx  f   dx' 


which  is  a  function  of  x  alone. 
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The  solution  of  equation  (3)  is  expressed  In  terras  of  modi- 
fied Bessel  functions  as 

l0(v)  -  l0(v  X) 

u   "   -f  '   qryT-^-  <*> 

where  Y  =  3rQ,  rQ  Is  the  radius  of  the  tube. 

Equation  (k)    will  define  the  velocity  profiles  In  the  en- 
trance region  If  the  relation  between  y  and  x  Is  known. 

When  3  =  0,  (that  is,  y  =  0),  the  acceleration  terms  are 
zero.   The  solution  then  assumes  the  form  of  a  fully  developed 
parabolic  profile.   When  3  — *■  oo,  (that  Is,  y  ~ "*  °°)  the  solution 
assumes  a  uniform  profile  at  the  entry.   Therefore,  the  solution 
provides  a  smooth  transition  from  the  entry  where  y  — >  oo  to  the 
fully  developed  section  where  y  =  0. 

For  determining  the  relation  between  y  and  x,  the  momentum 
integral  equation  is  used.   Integrating  the  boundary  layer  equa- 
tion in  cylindrical  coordinates  over  the  cross  section,  we  obtain 

2  A.   /Vrdr  -  K  r2  .  I  d£  r2  +  2  „  /r0(a£u   1  i*)rir       (O 
dx  f.      u  rar  -  gxr0   f  dx  r0  *  '"'      l„  2  *  r  3r'rar-   °' 
v  j  o   3r 

After  the  substitution  of  (*(•)  and  simplification,  we  obtain  the 
following  relations: 

1/2 
X  =  /    YdZ  (6) 

Z 


where 


x  ■  ^K  (7) 
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X  = *-. g  <•> 


EyI^y)  -  y' 


and 


^I0(y)I2(y)  -  [I0(y)  -  l]2  -  2I^(y) 

-6  =  p .        (9) 

2l2(y) 

It  can  be  shown  that  the  limits  for  variables  Z  and  X  ares 

Z  — »  j  <       X  — *  0     as     y  — *  oo 

2 
Z  — *■  ~  ?  t      ¥  — *■   oo     as     y  — *•  0  . 

Equation  (6)  provides  a  relation  between  x  and  y.  The 
Integration  can  be  carried  out  numerically  as  will  be  shown 
later  on. 

The  entrance  length  is  defined  as  the  distance  from  the 
entry  to  a  section  where  the  center  line  velocity  is  99  per  cent 
of  the  fully  developed  one;  that  Is, 

uQ  =  0.99  (u0)p>D> 


or 


Since 


_°-  =  o.99  -■  F•P• 


(Vf.d.  m    2 


hence 
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u0 
UQ  =  -^     =  0.99  x  2  =  1.98  . 


This  is  equivalent  to  saying  that  the  entrance  length  Is  the  dis- 
tance from  the  entry  to  a  section  where  the  dimenslonless  center 
line  velocity,  UQ,  Is  equal  to  1.98.   The  values  of  UQ  at  various 
itatned  from  flj.}  with 


sections  can  be  obtained  from  (k)    with  —  =  0|  that  Is, 


uQ     IQ(Y)  -  1 

uo  -  f  =  -\r^~  '  (10) 

Hence,  when  UQ  =  1.98,  the  corresponding  value  of  v  (and  there- 
fore the  value  of  x)  can  be  determined. 

In  calculating  the  pressure  drop  in  the  entrance  region, 
Langhaar  applied  the  mechanical  energy  balance  Instead  of  using 

the  Bernoulli  equation  because  in  the  downstream  part  of  the 

region 
entrance^where  the  velocity  approaches  the  fully  developed  one, 

the  Bernoulli  equation  ceases  to  be  applicable.   The  energy  loss 

due  to  viscous  dissipation  was  considered.   Equating  the  work  of 

pressure  force  and  frlctlonal  force  to  the  change  in  kinetic 

energy  gives  the  energy  equation 

x        2 
uA(p  -  p)   =  Ml   dx  /(!£)  dA  +  i  f  /(u3  -  53)dA       (11) 
0    A  dr       A        A 

where  A  is  the  cross-sectional  area  of  the  tube  and  M  is  the 
viscosity  of  the  fluid. 

Substituting  (k)    into  (11)  and  simplifying,  gives  the  total 
pressure  drop  between  a  supply  reservoir  and  the  section  x  in 
the  entrance  region: 
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P    = 


where 


P-P  Xv2,  1, 

«     1+2  /     *g  (IJ  -  I0I2)dX  +  2  /    (U^-DHdR       (12) 


|.fa2  o   i^ 


u   -   -a.  ,  a  »  x.. 

u  r0 

Equation  (12)  can  be  carried  out  numerically  and  will  be 
shown  in  the  later  part  of  this  chapter.   Hence  the  pressure 
drops  in  the  entrance  region  are  defined. 

The  detailed  derivation  of  the  equations  is  given  in  Li's 
report  [2], 
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3.2.  Computational  Procedures 

The  available  modified  Bessel  function  data  for  arguments 
from  0  to  85  [  3 ]  are  used  In  this  computation.   The  computation  - 
for  entrance  length,  velocity  distribution  and  pressure  distribu- 
tion will  be  given  separately  by  three  different  FOHTKAN  pro- 
grams. 

(1)   Entrance  Length 

The  relationship  between  the  length  in  the  entrance  region, 
X,  and  the  parameter,  Y,  Is  given  by  equation  (6).   The  function 
X(y)  is  represented  by  the  area  under  a  curve  (see  Fig.  3.1) 
represented  by  X(y).  which  is  a  function  of  Y  as  shown  in  equa- 
tion (0),  in  the  interval  of  Z(y),  (Z,  j) ,  which  is  also  a 
function  of  Y  as  shown  In  equation  (9).   By  assigning  a  set  of 
values  for  y,  a  set  of  the  corresponding  values  for  Z(y),  and 
Y(y)  can  be  obtained  from  equations  (9)  and  (8)  respectively. 
Then  the  corresponding  values  of  X(y)  are  obtained  by  numerical 
integration  of  equation  (6). 

The  set  of  values  assigned  to  y  are  chosen  from  100  to  0.1 
and  the  Intervals  between  each  y  is  shown  In  Table  3.1. 

The  values  of  Iq(y)  and  it(y)  for  Y  =  100>  95>  9°  are  ob- 
tained from  the  following  polynomial  approximate  relations  [6], 

1 
Y  e  yIq(y)  =  3.989W28  +  0.1328592t_1  ♦  .00225319t"2 


.00157565t"3  +  .009l628lt"4'  -  .02057706t"5 
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Table  3.1.   Selected  Interval  between  each  y 
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Y 

AY 

N* 

100  . 

5.0 

15 

30  : 

1.0 

15 

15  : 

0.5 

16 

7.0  : 

0.2 

20 

3.0  : 

0.1 

29 

0.1  -1 

N»:   Points  within  the  interval. 


3^ 

+  .02635537t"6  -  .Ol647633t" 7 
f   ,00392377t" 8  +  €  .  (13) 

where 

3.75  i  y  i  -  •       t  -  ^j 

)fc|  <  1.9  x  10"7  , 
and 

1 
Y  e"YI1(Y)   =   .39894228  -  .03988024-t"1  -  ,003620l8t"2 

■»■  .ooi638oit-3  -  .0l03i555t_i+  +  .02282967t"5 

-  .02895312t"6  +  .0178?654t"7 

-  .oo^2oo59t" 8  +  6  .  (ik) 

where 

3-75  £  y  £  «  .       t  =  y^j 

\€|  <  2.2  x  10"7  . 

In  the  computation,  the  values  of  I0(y)  and  I-i(y)  for  tne 
corresponding  values  of  y  are  given  as  the  input  data  In  the 
FOftTBAN  program  and  the  value  of  I2(y)  corresponding  to  each 
value  of  y  is  computed  by  the  following  relation  [  fcji 

X2(^  =  "  f  xi(y>  ♦  jo(y)  115) 

The  Integration  of  equation  (6)  proceeds  from  the  large 
value  of  y  to  the  small  one  and  the  trapezoidal  method  [  5  ]  is 
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applied  In  the  integration  (see  Fig.  3.1).   The  center  line 
velocity  corresponding  to  each  entrance  length  is  also  computed 
from  equation  (10). 

A  flow  chart  (Table  3.2)  will  be  followed  by  a  symbol  table 
(Table  3.3)  and  FORTRAN  program  (Table  J.k) .      The  whole  computa- 
tion is  carried  out  by  an  IBM  1620  computer  and  some  numerical 
results  are  tabulated  in  Table  3»5« 
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Table  3.2.     Flow  chart  of  Langharr's  Method 

—  entrance  length 

ISTAM) 


Zm-S,  Ym'O 

X'O 


Compute 
Yn  (   E$-  8  J 


I 


Coma*  te 
t>/   Trapezoidal 

r«'e 


Table   3.3.      Symbol  Table 
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FORTRAN  Symbol 

Variables 

BIO(N) 

W 

BI1(N) 

i^y) 

GAilA 

Y 

ZH 

Zm 

YM 

*m 

DELTG 

A  Y 

BI2N 

i2(y) 

ZN 

Zn 

YN 

*n 

DZ 

*  Z 

X 

X 

uo 

°0 
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TABLfc  3.4.   FORTRAN  PROGRAM  FOR  LANGHAAr's  METHOD 
---ENTRANCE  LENGTH 


c 

LANGHAAR.S  METHOD  ENTRANCE  LENGTH 

n  1 

FORMATf 1E12.4) 

in 

F0PMAT(6F12.4) 

1  l 

FORMAT! 12. F5.1 .  2E11.4,  4E10, 

DIMFNSION  BIO(IOO).  BIK100) 

READ  lll.(BIO(N) .N=l .3) 

READ  lU.tBIOINI »N=4>95) 

READ  111 • <  BI 1 (N) »N=1 .3) 

READ   lO.lBIKN)  »N  =  4.95> 

N  =  0 

GAMA=100. 

ZM=.5 

YM=.0 

x  =  .o 

DELTG=.0 

.3) 

? 

GAMA=GAMA-DFLTG 

M  =  N  + 1 

BI2N=-?.»BI1 (N)/GAMA+B!0(N) 

ZN=2.*BI0(N)/BI2N-((B!0(N)-1. 

,  )/BI2N)»( (BlO(N)-l 

ZN=ZN-(BI1(N)/BI2N)»(BI1(N)/BI2N) 

YN=BI2N/(2.*GAMA*BI1 ( N ) -GAMA»GAMA ) 

DZ=ZM-ZN 

20 

X=X+(YM+YN)»DZ/2. 
U0=BIC(N)/BI2N-1./BI2N 

PUNCH  11,  N.GAMA.X.UO.BIOIN) , 

.BIKN1.ZN.YN 

ZM  =  ZN 

YM  =  YN 

IFIGAMA-30. 131.31. 3 

3 

DELTG=5. 
GO  TO  2 

31 

IFIGAMA-15. 141.41.4 

4 

DELTG=1. 
GO  TO  2 

&] 

1FIGAMA-7. 151.51.5 

5 

TFLTG=.5 
GO  TO  2 

Ri 

IF(GAMA-3.)61.61.6 

6 

DFLT6».2 

GO  TO  2 

61 

IF1GAMA-.1 171.71.7 

7 

DFLTG  =  .'. 
GO  TO  2 

71 

STOP 
END 

i  )/BI2N)/2. 
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TABLE  3.5.   NUMERICAL  RESULTS  CF  LANGHAAR's  METHOD ENTRANCE  LENGTH 

N  GAMA      X  UO        RIO(N)     BIl(N)       ZN        YN 

1100.0  0.2772E-04  0.1020E+01  0.107E+43  0.107E+43  0.489  0.492E-02 

2  95.0  0.2570E-04  0.1021E+01  0.742E+40  0.738E+40  0.489  0.518E-02 

3  90.0  0.2960t-QA  0.1023E+01  0.514E+38  0.511E+38  0.488  0.546E-02 

4  85.0  0.3463E-04  0.1U24E+01  0.356E+36  0.354E+36  0.488  0.578E-02 

5  80. 0  0.4213E-04  0.1025E+01  0.248E+34  0.246E+34  0.486  0.613E-02 

6  75.0  0.4161E-04  0.1027E+01  0.172E+32  0.171E+32  0.486  0.654E-02 

7  70.0  0.4834E-04  0.1029E+01  0.120E+30  0.119E+30  0.485  0.699E-02 

8  65.0  0.6229E-04  0.1031E+01  0.840E+27  0.834E+27  0.483  0.752E-02 

9  60.0  0.7517E-04  0.1034E+01  0.r89E+25  0.585E+25  0.482  0.813E-02 
in  55.0  0.8892E-04  0.1037E+01  0.*15E+23  0.411E+23  0.480  0.884E-02 

11  50.0  0.1061E-03  0.1041E+01  0.293E+21  0.290E+21  0.478  0.970E-02 

12  45.0  0.1410E-03  0.1046E+01  0.208E+19  0.206E+19  0.475  0.107E-01 

13  40.0  0.1900E-03  0.1052E+01  0.149E+17  0.147E+17  0.471  0.120E-01 
U  35.0  0.2473E-03  0.1060E+01  0.107E+15  0.106E+15  0.466  0.137F-01 

15  30.0  0.3253E-03  0.1070E+01  0.782E+12  0.769E+12  0.461  0.158E-01 

16  29.0  0.3522E-03  0.1073E+01  0.293E+12  0.287E+12  0.459  0.164E-01 

17  28.0  C.4088E-03  0.1075E+01  0.11CE+12  0.108E+12  0.456  0.169E-01 

18  27.0  0.4110E-03  0.1078E+01  0.410E+11  0.403E+11  0.456  0.175E-01 

19  26.0  C.4471E-03  0.1082E+01  0.154E+H  0.151E+11  0.454  0.181E-01 

20  25.0  0.4892E-03  0.1085E+01  0.578E+10  0.566E+10  0.451  0.188E-01 

21  24.0  0.5317E-03  0.1089E+01  0.217E+10  0.212E+10  0.449  n,195j>01 

22  23.0  0.5942E-03  0.1093E+01  0.815E+09  0.797E+09  0.446  0.203E-01 

23  22.0  0.6475E-03  0.1097E+01  0.307E+09  0.300E+09  0.443  0.212E-01 

24  21.0  0.7210E-03  0.H02E+01  0.116E+09  0.113E+09  0.440  0.221E-01 

25  20.0  0.8051E-03  0.1108E+01  0.436E+08  0.425E+08  0.436  0.232E-01 

26  19.0  0.9111E-03  0.1114E+01  0.165E+08  0.160E+08  0.432  0.243E-01 

27  18.0  0.1039E-02  0.1121E+01  0.622E+07  0.604E+07  0.427  0.255E-01 

28  17.0  0.1192E-02  0.1129E+01  0.236E+07  0.229E+07  0.421  0.269E-01 

29  16.0  0.1366E-02  0.H38F+01  0.893E+06  0.865E+06  0.415  0.284E-01 

30  15.0  0.1585E-02  0.1148E+01  0.340E+06  0.328E+06  0.407  0.301E-01 

31  14.5  0.1694E-02  0.1I53E+01  0.210E+06  0.202E+06  0.403  0.310E-01 

32  14.0  0.1865E-02  0.H60E+01  0.129E+06  0.125E+06  0.398  0.320E-01 

33  13.5  0.2014E-02  0.1166E+01  0.800E+05  0.769E+05  0.393  0.330E-01 

34  13.0  0.2217E-02  0.1173E+01  0.494E+05  0.475E+05  0.387  0.341E-01 

35  12.5  0.2435E-02  0.1181E+01  0.306E+05  0.294E+05  0.381  0.353E-01 

36  12.0  0.2669E-02  0.1190E+01  0.190E+05  0.181E+05  0.375  0.366E-01 

37  11.5  0.2905E-02  0.1199E+01  0.118E+05  0.112E+05  0.368  0.380E-01 
3fl  11.0  0.3324E-02  0.1210E+01  0.729E+04  0.695E+04  0.357  -  0.394E-01 

39  10.5  0.3718E-02  0.1221E+01  0.453E+04  0.431E+04  0.348  0.410E-01 

40  10.0  0.4165E-02  0.1234E+01  0.282E+04  0.267E+04  0.337  0.428E-01 


<+0 


N  GAMA      X           UO        BIO(N)     BIKN)       ZN        YN 

41  9.5  0.4737E-02  0.1248E+01  0.175E+04  0.166E+04  0.324  0.447E-01 

4?  9.0  0.5363E-02  0.1264E+01  0.109E+04  0.103E+04  0.310  0.468E-01 

43  8.5  0.6208F-02  0.1282E+01  0.683E+03  0.642E+03  0.293  0.491E-01 

44  8.0  0.7214E-02  C.1302E+01  0.428E+0?  0.400E+03  0.273  0.517E-01 

45  7.5  0.8425E-02  0.1325E+01  0.268E+03  0.250E+03  0.250  0.547E-01 

46  7.0  0.9903E-02  0.1351E+01  0.169E+03  0.156E+03  0.224  0.581E-01 

47  6.8  0.1079E-01  0.1363E+01  0.140E+03  O.129E+03  0.209  0.595E-01 

48  6.6  0.1153E-01  0.1375E+01  0.117E+03  0.107E+03  0.196  0.612J-01 

49  6.4  0.1233E-01  0.13B8E+01  0.970E+02  0.890E+02  0.183  0.629E-01 

50  6.2  0.1327E-01  0.1401E+01  0.807E+02  0.739E+02  0.169  0.648E-01 

51  6.0  0.1431E-01  0.1416E+01  0.672E+02  0.613E+C2  0.153  0.668E-01 
5?  5.8  0.1546E-01  0.1431E+01  0.560E+02  0.510E+C2  0.136  0.690E-01 

53  5.6  0.1672E-01  0.1446E+01  0.467E+02  0.423E+02  0.118  0.714E-01 

54  5.4  0.1809E-01  0.1463E+01  0.390E+02  O.352E+02  0.992E-01  0.741E-01 

55  5.2  0.1967E-01  0.1481E+01  0.326E+02  0.293E+02  0.784E-01  0.770E-01 

56  5.0  0.2149E-01  0.1499E+01  0.272E+02  0.243E+02  0.552E-01  0.801E-01 

57  4.8  0.2334E-01  0.1518E+01  0.228E+02  0.203E+02  0.326E-01  0.838E-01 

58  4.6  0.2546E-01  0.1138E+01  0.191E+02  0.169E+02  0.793E-02  0.878E-01 

59  4.4  0.2802E-01  0.1i60E+01  0. 160E+02  0. 141 E+02-0. 206E-0 1  0.923E-01 

60  4.2  0.3086E-01  0.1582E+01  0.134E+02  0. 117E+02-0. 505E-01  0.974E-01 

61  4.0  0.3354E-01  0.1604E+0]  0.113E+02  0.976E+01 -0. 77 1E-01  0.103 
6?  3.8  0.3672C-01  0.1628E+01  0.952E+01  0. 814E+01-0. 107  0.110 

63  3.6  0.4O52E-01  C.1652E+01  0.803E+01  0.679E+01-0 . 140  0.118 

64  3.4  0.4467E-01  0.1677E+01  0.679E+01  0. 567E+01-0. 174  0.128 

65  3.2  0.4944E-01  0.1703E+01  0.575E+01  0.473E+01-0. 2 10  0.139 

66  3.0  0.5458E-01  0.1728E+01  0.488E+01  0.395E+01-0. 245  0.153 

67  2.9  0.5787E-01  0.1742E+01  0.450E+01  0. 361 E+01-0. 266  0.160 

68  2.8  0.6098E-01  0.1755E+01  0.M6E  +  01  0. 330E+01-0. 285  0.169 

69  2.7  0.6388E-01  0.1767E+01  0.i84E+01  0. 302E+01 -0.30 1  0.179 
7n  2.6  0.6753E-01  0.1781E+01  0.355E+01  0. 276E+01-0 . 321  0.190 

71  2.5  C.7175E-01  0.1794E+01  0.329E+01  0.252E+01 -0. 343  0.201 

7?  2.4  0.7562E-01  0.1807E+01  0.305E+01  0. 230E+01-0. 36 1  0.215 

73  2.3  0.7969E-01  0.1820E+01  0.283E+01  0.210E+01-0. 380  0.231 

74  2.2  0.8446E-01  C.1832E+01  0.263E+01  0. 191 E+01-0. 400  0.248 
7?  2.1  0.8823E-01  0.1844E+01  0.245E+01  0. 1 75E+01-0.41 4  0.269 

76  2.0  0.9530E-01  0.1858E+01  0.228E+0]  0. 159E+01-0.440  0.291 

77  1.9  0.9776E-01  C.1868E+01  0.213E+01  0. 145E+01-0.448  0.319 

78  1.8  0.1029  0.1880E+01  0.199E+0]  0. 132E+01-0.463  0.351 

79  1.7  0.1092  0.1891E+01  0.186E+01  0. 120E+01-0.480  0.388 

80  1.6  0.1249  0.1905E+01  0.175E+01  0. 109E+01 -0. 518  0.432 


N  GAMA      X           UO  BIO(N)  BIKNI  ZN        YN 

0.1914E+01  0.165E  +  01  0.982  -0.520 

0.1926E+01  0.155E+01  0.886  -0.560 

0.1935E  +  01  0.147E+01  0.797  -0.571 

0.1942E+01  0.139E+01  0.715  -0.557 

■0.1953E+01  0.133E+01  0.638  -0.607 

0.1">62E  +  01  0.127E  +  01  0.565  -0.625 

87  .9  0.1653      0.1966E+01  0.121E+01  0.497  -0.59* 

88  .8  0.1868E-01  0.1971E+01  0.117E+01  0.433  -0.493 

89  .7  0.6694      0.1986E+01  0.113E+01  0.372  -0.847 

90  .6  0.2484      0.1986E+01  0.109E+01  0.314  -0.674 

91  .5  0.3954E+01  0.2U06E+01  0.106E+01  0.258  -0.177E+01  0.397E+01 

92  .4  0.5358E+01  0.2000E+01  0.104E+01  0.204  -0.204E+01  0.625E+01 

93  .3-0.5119E+02  0.1971E+01  0.102E+01  0.152  0.435E+01  0.114E+02 

94  .2  0.6476E+02  0.2000E+01  0.101E+01  0.101  -0.201E+01  0.250E+02 

95  .1-0.2983E+05  0.1667E+01  0.100E+01  0.501E-01  0.340E+03  0.150E  03 
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0, 

i282E+01 

k2 
(2)  Velocity  Distribution 
Equation  (k) 

Vy)  -  i0(v  *) 

u  "  "t  =  W-^  c*) 


represents  the  velocity  distribution  at  any  value  of  y  In  the 


entrance  region.   As  shown  in  the  previous  article  y  represents 

a  corresponding  value  of  x,  the  distance  from  the  entry.   The 

velocity  distribution  Is  assumed  to  be  symmetric  with  respect  to 

the  center  line.   The  values  of  the  velocity  at  eleven  points 

across  the  half  circular  tube,  (-f-  =  1.0,  0.9 0.1,  0)  at  a 

r0 
given  value  of  y  are  computed  from  equation  (k) .   Table  3.6  shows 

the  values  of  y  and  y  tt-.  where  the  velocity  is  computed. 
r0 
A  flow  chart  (Table  3.7)  for  this  computation  is  followed  by 

a  symbol  table  (Table  3.8)  and  POHTRAN  program  (Table  3.9).   The 

computation  is  also  carried  out  by  an  IBM  1620  computer  and  some 

numerical  results  are  tabulated  in  Table  3.10  and  the  velocity 

distribution  for  each  location  in  the  entrance  region  is  plotted 

as  shown  In  Pig.  3.2. 
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Table  3.7      Flow  chart  of  Langharr's  Method 

—  velocity  distribution 
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Table   3.8.      Symbol  Table 


^5 


FORTRAN  Symbol 

Variables 

BIO(N) 

Vy) 

BI1(N) 

^(y) 

BIOO(J) 

h&0  v) 

GAMA 

Y 

DELTG 

A  Y 

R 

r 
r0 

BI2N 

i2(y) 

U 

u 
u 

46 


TABLE  3.9.   FORTRAN  PROGRAM  FOR  LANGHAAr's  METHOD 
—  VELOCITY  DISTRIBUTION 


C   LANGHAmR.S  METHOD VELOCITY  DISTRIBUTION 

10  F0RMATI6E12.4) 

11  F0RMAT(2I4,2F6.1»4E12.4) 

12  F0RMAT(4X,I4,6X,F6.1,2E12.4) 
DIMENSION  BI0(12),BI1(12)»BI00(120) 

100  READ  10,  (BIO(N).  N=1.10) 

101  READ  10,  (BIKN),  N=l,10) 

102  READ  10,  (BIOO(J),  J=l,110) 
N  =  0 

J  =  0 

GAMA=40. 

DELTG=.0 

2  GAMA=GAMA-DELTG 
R  =  l. 

N  =  N  +  1 
J  =  J  +  1 

BI2N=-2.*BI1(N)/GAMA+BI0(N) 
20  U=(BI0(N)-BI0O( J) 1/BI2N 
IFIR-1.)  30,3,3 

3  PUNCH  11,  N,J,GAMA,R,U,BI00( J) ,BI0(N) ,BI1 (N) 
GO  TO  31 

30  PUNCH  12,  J,R,U,BI00(J) 
IFIR-O.)  32,32.31 

31  J=J+1 
R=R-.l 
GO  TO  20 

32  IFIGAMA-20.)  41,41,4 

4  DELTG=20. 
GO  TO  2 

41  IFtGAMA-15.)  5,5,42 

42  DELTG=5. 
GO  TO  2 

5  IFiGAMA-8.)  6,6,51 
51  DE*."G  =  7. 

GO  re   2 

6  IFiGAMA-6.)  7,7,61 
61  DELTG=2. 

GO  TO  2 

7  IFIGAMA-1.)  8,8,71 
71  DELTG=1. 

GO  TO  2 

8  STOP 
END 


<*7 


TABLF    3.10.       NUMERICAL    RFSULTS    OF    LANGHAAR'S    METHOD 
-—VELOCITY    DISTRIBUTION 


J   GAMA 

R 

U 

BIOO(J) 

BIO(N) 

Br i(n> 

1   40.0 

1.0 

0.0000 

0.1489E+17 

0.1489E+17 

0.1471E+17 

2 

.9 

0.1032E+01 

0.2877E+15 

3 

.6 

0.1052E+01 

0.5591E+13 

4 

.7 

0.1052F+01 

0.1095E+12 

5 

.6 

0.1052E+01 

0.2169E+10 

6 

.5 

0.1052F+01 

0.4356E+08 

• 

7 

.4 

0.1052E+01 

0.8934E+06 

8 

.3 

0.1052E+01 

0.1895E+05 

9 

.2 

0.1052E+01 

0.4276E+03 

10 

.1 

1.1052E+01 

0.1130E+02 

11 

0.0 

J.1052E+01 

0.1000E+01 

12   20.0 

1.0 

0.0000 

0.4356E+08 

0.4356E+08 

0.4245E+08 

13 

.9 

0.9498 

0.6218E+07 

14 

.8 

1.1085E+01 

0.8934E+06 

15 

.7 

0.1105E+01 

0.1294E+06 

16 

.6 

0.1107E+01 

0.1895E+05 

1  7 

.5 

0.1108E+01 

0.2816E+04 

18 

.4 

0.1108E+01 

0.4276E+03 

19 

.3 

0.1108F+01 

0.6723E+02 

20 

.2 

0.1108E+01 

0.1130E+02 

2  1 

.1 

0.1108E+01 

0.2280E+01 

22 

0.0 

0.1108E+01 

0.1000E+01 

23   15.0 

1.0 

0.0000 

0.3396E+06 

0.3396E+06 

0.3281E+06 

24 

.9 

0.8776 

0.7997E+05 

?5 

.8 

0.1084E+01 

0.1895E+05 

26 

.7 

0.1 133E+01 

0.4527E+04 

27 

.6 

0.1144E+01 

0. 1094E+04 

28 

.  5 

0.1147E+01 

0.2682E+03 

29 

.4 

0.1148E+01 

0.6723E+02 

30 

.3 

0.1148E+01 

0.1748E+02 

31 

.2 

0.1148E+01 

0.4881E+01 

32 

.1 

0.1148E+01 

0.1647E+01 

33 

0.0 

0.1148E+01 

0.1000E+01 

34    8.0 

1.0 

0.0000 

0.4276E+03 

0.4276E+03 

0.3999E+03 

35 

.9 

0.6858 

0.2029E+03 

36 

.8 

0.1009F+01 

0.9696E+02 

37 

.7 

0.U62E+01 

0.4674E+02 

38 

.6 

0.1236E+01 

0.2279E+02 

<+8 


J   GAf 

4A      R 

u 

BIOOI J) 

BIO(N) 

B  1 1  (  N  ) 

30 

.5 

0.1271E+01 

0.U30E+O2 

40 

.4 

0.1288F+01 

0.5747E+01 

41 

.3 

0.1296E+01 

0.3049E+01 

4? 

.2 

0.1300E+01 

0.1750E+01 

4' 

.1 

0.1302E+01 

0.1167E+01 

4/. 

O.'J 

0.1302E+01 

0.1000E+01 

45     6. 

.0    1.0 

0.0000 

0.6723t+u2 

0.6723E+02 

0.6134E+02 

46 

.9 

0.6032 

0.3901E+02 

47 

.8 

0.9499 

0.2279E+02 

4  8 

.7 

1.1150E+01 

0.1344E+02 

40 

.6 

J.1265E+01 

0.8028E+01 

SO 

.5 

0.1333F+01 

0.4881E+01 

51 

.4 

0.1372E+01 

0.3049E+01 

52 

.3 

O.1395E+01 

0.1990E+01 

53 

.2 

0.1407E+01 

0.1394E+01 

S4 

.1 

0.1414F+01 

0.1092E+01 

55 

0.0 

0.1416E+01 

0.1000E+01 

56    5. 

.0    1.0 

0.0000 

0.2724E+02 

0.2724E+02 

0.2434E+02 

57 

.9 

0.5576 

0.1748E+02 

5R 

.8 

0.9106 

0.1130E+02 

59 

.7 

0.1135E+01 

0.7378E+01 

60 

.6 

0.1277E+01 

0.4881E+01 

61 

.5 

0.1368E+01 

0.3290E+01 

62 

.4 

0.1426E+01 

0.2280E+01 

63 

.3 

0.1462E+01 

0.1647E+01 

6  4 

•  2 

0.1484E+01 

0.1266E+01 

6  5 

.  1 

0.1495E+01 

0.1063E+01 

66 

0.0 

0.1499E+01 

0.1000E+01 

67    4. 

0    1.0 

0.0000 

0.1130E+O2 

0.1130E+02 

0.9759E+01 

68 

.9 

0.5096 

0.8028E+01 

60 

.8 

0.8649 

0.5747E+01 

70 

.7 

0.1113E+01 

0.4157E+01 

71 

.6 

0.1285F+01 

0.3049E+01 

72 

.5 

0.1405F+01 

0.2280E+01 

73 

.4 

0.1487E+01 

0.1750E+01 

74 

.  3 

0.1543E+01 

0.1394E+01 

75 

.2 

0.1578E+01 

0.1167E+01 

76 

.1 

0.1598E+01 

0.1040E+01 

77 

0.0 

0.1604E+01 

0.1000E+01 

<*9 


N    J   GA» 

1A            R 

u 

BI00(J) 

BIO(N) 

RUIN) 

8   78    3. 

.0    1.0 

0.0000 

0.4881E+01 

0.4881E+01 

0.3953E+01 

79 

.9 

0.4627 

0.3842E+01 

80 

.8 

0.8158 

0.3049E+01 

3  1 

.7 

0.1084E+01 

0.2446E+01 

82 

.6 

0.1287E+01 

0.1990E+01 

83 

.5 

0.1440F+01 

0.1647E+01 

«4 

.4 

0.1553F+01 

0.1394E+01 

85 

.3 

0.1633F+01 

0.1213E+01 

86 

.2 

0.1687F+01 

0.1092E+01 

37 

.1 

0.1718F+01 

0.1023E+01 

RB 

0.0 

0.1728E+01 

0.1000E+01 

9   89    2. 

.0    1.0 

0.0000 

0.2280E+01 

0.2280E+01 

0.1591E+01 

90 

.9 

0.4209 

0.1990E+01 

91 

.8 

0.7692 

0.1750E+01 

92 

.7 

0.1055E+01 

0.1553E+01 

93 

.6 

0.1286E+01 

0.1394E+01 

94 

.? 

0.1472F+01 

0.1266E+01 

95 

.4 

0.1615E+01 

0.1167E+01 

96 

.3 

0.1724E+01 

0.1092E+01 

97 

.2 

0.1800F+01 

0.1040E+01 

9  P. 

.1 

0.1843E+01 

0.1010E+01 

99 

0.0 

0.1858E+01 

0.1000E+01 

10  100    1. 

0    1.0 

0.0000 

0.1266E+01 

0.1266E+01 

0.5652 

101 

.9 

0.3909 

0.1213E+01 

102 

.8 

0.7301 

0.1167E+01 

103 

.7 

0.1032E+01 

0.1126E+01 

104 

•  6 

0.1283E+01 

0.1092E+01 

105 

.5 

0.1497F+01 

0.1063E+01 

106 

.4 

0.1667F+01 

0.1040F+H1 

107 

.3 

0.1792E+01 

0.1023E+01 

. 

108 

.2 

1.1888E+01 

0.1010E+01 

109 

.  1 

U.1940E+01 

0.1003E+01 

110 

0.0 

0.1962E+01 

0.1000E+01 
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51 
(3)  Pressure  Distribution 

Equation  (12)  represents  the  total  pressure  drop  between  the 
supply  reservoir  and  section  x  in  the  entrance  region. 
Let 

*(y)   =  ^  [ij  -  I0I2]  (16) 


then  we  have 
X 


;    Tf  IX1  "  I0I2ldX  =   '   *(y>«  •  (1?) 


The  evaluation  of  this  integral  is  difficult  due  to  the  fact  that 
0(y)  — *  w  as  x  — »  0  (y  — *•   oo)  .   Therefore  we  make  use  of  a 
limiting  relation  (see  p.  100  In  [2])  and  obtain  a  result  as 
followsi 


\  I     (y"3  +  ...)4y  =  -1* 

Y  ^Y2 


It  follows  that 


lim  Y2X  =  | 


It  is  pointed  out  by  Langhaar  that  setting  y2X  =  i  In  the  range 
Y  >  100  is  a  good  approximation.   With  this  approximation  the 
integration  of  equation  (1?)  becomes  easy  and  is  carried  out  as 
follows: 


52 

X  .00002772  X 

/  0(y)dX  =  /         0(y)dX  +  /        0(y)dX 
0  0  .00002772 

where  X  =  .00002772  corresponds  to  y  =  100.  Furthermore,  we  can 
show  (see  [2],  p.  10^) 

I2   I 

I2    h  V        y2  +  y3         ' 

Making  the  approximation  that  y  X  =  £  for  y  2  100,  we  obtain 

I2   I 
*(Y)   =  Y2  (4  -  A 
l|  h 


~  Y  ♦ 


3  +  £ 


1        1 

1     2  ? 

~  ±  X   "  *  3  +  lOX4 


and 


.00002772  .00002772  ,   -  i         i 

/         0(y)dX  =  /         (±X   2  +  3  +  10X2)dX 

0  0  c 

=  .00535  (18) 


Therefore, 

X  X 

/  *(y)dX  =  .00535  ♦  /        0(y)dX 

0  .00002772 

so  that  the  numerical  evaluation  of  equation  (12),  the  pressure 
distribution,  proceeds  as  follows: 
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X 

(i)   Evaluating  /  0(y)dX:   The  integral  here  represents 

0 
the  area  under  the  curve  In  Fig.  3.3,  part  (b).   The 

range  of  X  is  0  to  »,  corresponding  to  the  value  of  y 
from  °o  to  0,  and  the  range  of  0(y)  is  0  to  »,  corre- 
sponding to  the  value  of  y  from  0  to  «>.   Since  the  area 
under  the  curve  corresponding  to  y  from  »  to  100  has 
been  evaluated  as  shown  in  equation  (18),  we  will  put 
Y  =  100  as  the  Initial  value  and  assign  a  set  of  y  as 
indicated  in  Table  3.1.   Then  the  X  value  for  each  y 
can  be  computed  by  using  equation  (6)  (also  see  Fig. 
3. 3 i  part  (a)),  and  the  value  of  0(y)  can  be  computed 
from  equation  (16).   The  trapezoidal  method  is  also 
applied  in  this  numerical  integration. 

1   o 
(ii)  Evaluating  /  (UJ-l)RdR:   For  a  given  value  of  y, 

0 
this  integral  represents  the  area  under  the  curve  in 

Fig.  J.k.      The  range  of  R  is  0  to  1.   We  will  divide  the 

range  of  B  into  20  equi -distant  intervals,  so  that  the 

area  is  divided  into  20  longitudinal  strips  and  the 

total  area  is  the  sum  of  the  areas  of  these  strips. 

The  area  of  each  strip  is  evaluated  by  applying  the 

trapezoidal  method  and  the  initial  value  is  chosen  at 

the  origin.   The  value  of  (iP-l)R  at  a  given  value  of 

y  can  be  computed  by  using  equation  (4)  for  different 

values  of  H. 
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(lil)  As  soon  as  the  above  two  Integrals  have  been 
evaluated,  the  total  pressure  drop  between  the  supply 
reservoir  and  section  X  in  the  entrance  region  corre- 
sponding to  a  certain  value  of  y  °an  be  computed  from 
equation  (12). 

A  flow  chart  (Table  3.10)  will  be  followed  by  a 
symbol  table  (Table  3.11)  and  FORTRAN  program  (Table 
3.12).   The  whole  computation  is  carried  out  by  using 
an  IBM  1620  digital  computer.   It  can  be  seen  that  in 
computing  the  pressure  drop,  the  entrance  length  and 
the  velocity  distribution  are  computed  at  the  same 
time.   The  statements  from  *+30  to  k?   In  this  program 
are  inserted  in  order  to  punch  out  the  different 
velocity  values  at  different  values  of  R  when  y  Is 
equal  to  8  and  7.5.   The  results  are  shown  In  Table 
3.13  which  can  be  compared  with  Table  3.5  and  Table 
3.10.   The  pressure  distribution  can  be  plotted  from 
the  results  and  Is  shown  in  Fig.  3.5. 
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Fig. 3.4.        (U3-l)RdR    represents  the    area    under 
Jo 

the  curve   at  a   given    value  of    T  . 
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Table  3.10.    Flow  chart  of  Langharr's  Method 

—  pressure  distribution 


a-- 1 


>r---2 
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('  (t*od      \ 
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Zf'48874 
Y«=.  004.9  D39 

<p(rL*  110.53+ 

Xm 

<f(f)d  x-;oo5iS 


i 


'—  r '  r~*r 


Compute 

X* 
<P(f)dx 

o 


f 


R-0 
f(R)m«0 


J 


^jT=.l 


-fsropr- "(end 


Com  pufe 

I.(  ru) 

17,    f(R)n 

PR 


f(ft)m'  f(R)n 


|  R'ff+.o5| 


<ro 


Zm:  Zr> 
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Table   3.11.      Symbol  Table 
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FORTRAN  Symbol 

Variable 

FORTRAN  Symbol 

Variable 

BIO(N) 

Vy) 

FGN 

x„ 

BI1(N) 

ix(y) 

DPX 

Xm 

GAMA 

Y 

R 

r/r0 

DELTG 

A  Y 

PR 

1   1 

/  (IK-l)RdR 
0 

ZM 

Zm 

FRM 

(U3-l)Rm 
in 

YM 

*m 

H 

Y-r/r0 

XM 

Xm 

T 

(y  t^)/3.75 
r0 

FGM 

E 

e  =  2.71828 

PX 

;  m^(Y)dx 

0 

BIOH 

I0(Y-r/r0) 

BI2N 

I2<Y> 

U 

u/u 

ZN 

Zn 

FRN 

(iP-l)^ 

IN 

Yn 

P 

(p„-p)/f/'s2 

DZ 

A  Z 

DX 

A  X 

XN 

Xn 

L 

6o 


TABLE  3.  12.   FORTRAN  PROGRAM  FOR  LANGHAAR'S  METHOD 
—-PRESSURE  DISTRIBUTION 

C   LANGHAAR.S  METHOD PRESSURE  DISTRIBUTION 

111  F0RMATI2E12.4) 

1  FCRMATI6E12.4) 

10  FORMAT(F5.1.E12.5.F7.3,E12.5) 

11  FeRMAT(F5.1.E12.5»F7.3.12Xi3E12.S> 

12  FORMAT! 5X.12X.F7.3.E12. 5. 36X1 
DIMENSION  BIO!  100)  .BIK100) 
RFAD  lll.(BIOtN) .  N=l,2) 
RFAD  Ii(BIOIN) t  N=3»94) 

RFAD  till  (BI1 <N) ,  N=1.2) 

RFAD  1  .  (BIKN)  .  N=3»94) 

N  =  0 

GAMA=95. 

DFLTG=.0 

ZM=. 488740 

YM=. 00492339 

XM=.277177E-04 

FGM=110.534 

PX=. 00535 

2  GAMA=GAMA-DELTG 
N  =  N+1 

BI?N=-7.*BI1 (Nl/GAMA+BIOIN) 

ZN=2.*BI0(N)/BI2N-< (BIO(N)-l. >/BI2N)#( (BIO(N)-l. ) /BI2N1/2. 

ZN=ZN-!BI1 <N>/BI2N)»(BI1(N)/BI2N> 

YN  =  BI2N/(2.*GAMA*BI1  (  N  )  -GAI'A#GAMA  ) 

DZ=ZM-ZN 

DX= (YM+YN)*DZ/2. 

XN=DX+XM 

3  FGN=GAMA*GAMA*( < BI 1 ( N > /BI 2N > *BI 1 ( N > /BI 2N  -BI 0 ( N ) /B I 2N ) 
DPX=(FGN+FGM)*DX/2. 

PX=PX+DPX 
31  R*.0 
PR  =  .0 
FRM=.0 

40  R=R+.05 

4  H=GAMA*R 
T=H/3«75 
IFIH-3.75)  41,41.4? 

41  BIOH=1.+3.5156229*T»T+3.08  99424*T**4.+1.206  74  92*T**6. 
BI0H=BI0H+.2659732*T*»8.+.O36O768*T*«lO.+.OO45813*T*»12. 
GO  TO  43 

42  E=2. 71828 

BI0H=.398942+.0132859*T»*I-1. ) +. 0022531 *T«* ( -2. ) -.001 5  756«T*» ( -3 . ) 
BI0H=BI0H+.00916281*T*«(-4. ) -.0205  770*T»« < -5. ) +.0263554*T»* ( -6. ) 
BI0H=(BIOH-.0164763*T**(-7. ) +. 0039237*T»* ( -8. ) ) *E**H/H*» I .5 ) 

43  U=(BI0(N) /RI2N1-BI0H/BI2N 
430  IFIGAMA-8. )44»44»431 
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44  1F(GADA-7.5>431,45,45 

45  IF(R-. 05)46. 46.47 

46  PUNCH  10,  GAMA.XN.R.U 
GO  TC  431 

47  PUNCH  12.  R.U 
431  FRN=(U*U*U-1. )»R 

PR=PR+(FRN+FRM)»l./40. 
FRM=FRN 

IFIR-1. 140,471.471 
471  P=1.+2.»(PX+PR) 

PUNCH  11,  GAMA.XN.R.PX.PR.P 

2M  =  ZN 

YM  =  YN 

XM  =  XN 

FGM=FGN 

IFIGAMA-30. 151,51,5 

5  DFLTG=5. 
GO  TC  2 

51  !F(GAMA-15.)61,61,6 

6  DFLTG=1. 
GO  TO  2 

61  !F(GAMA-7.)71,71,7 

7  DELTG=.5 
GO  TO  2 

71  IFIGAMA-3. 181.81,8 

8  DELTG=.2 
GO  TC  2 

81  IF(GAMA-.1)91,91 ,9 

9  DFLTG=.l 
GO  TC  2 

91  STOP 
FND 


TABLE  3.  13. 
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NUMERICAL  RESULTS  OF  LANGHAAR's  METHOD 
PRESSURE  DISTRIBUTION 


RAMA 


XN 


PX 


PR 


95.0  0.25698E-04  1.000 

90.0  0.29600E-04  1.000 

85.0  0.34632E-04  1.000 

80.0  0.42129E-04  1.000 

75.0  0.41607E-04  1.000 

70.0  C.48340E-04  1.000 

65.0  0.62291E-04  1.000 

60.0  0.75165E-04  1.000 

55.0  0.88917E-04  1.000 

50.0  0.10613E-03  1.000 

45.0  0.14100E-03  1.000 

40.0  0.18997E-03  1.000 

35.0  0.24726E-03  1.000 

30.0  0.32532E-03  1.0">0 

29.0  0.35222E-03  I.OjO 

28.0  0.40884E-03  1.000 

27.0  0.4110UE-03  1.000 

26.0  0.44713E-C3  1.000 

25.0  0.48918E-03  1.000 

24.0  0.5317OE-03  1.000 

23.0  0.59418F-03  1.000 

22.0  0.64747E-03  1.000 

21.0  0.72095E-03  1.100 

20.0  0.80507E-03  l.JOO 

19.0  0.91111E-03  1.000 

18.0  0.10386E-02  1.000 

17.0  0.11916E-02  1.000 

16.0  0.13659E-02  1.000 

15.0  0.15845E-02  1.000 

14.5  0.16941E-02  1.000 

14.0  0.18650E-02  1.000 

13.5  0.20142E-02  1.000 

1 3.0  0.22174E-02  1.000 

12.5  0.24347E-02  1.000 

12.0  0.26687E-02  1.000 

11.5  0.29051E-02  1.000 

11.0  0.33235E-02  1.000 

10.5  0.37175E-02  1.000 

10.0  0.41649E-02  1.000 

9.5  0.47368E-02  1.000 

9.0  0.53627E-02  1.000 

8.5  0.62083E-02  1.000 


0.51415E-02  0.47856E-02  0.10199E+01 

0.55078E-02  0.61964E-02  0.10234E+01 

0.59607E-02  0.77159E-02  0.10274E+01 

0.66132E-02  0.93582E-02  0.10319E+01 

0.65713E-02  0.11102E-01  0.10353E+01 

0.70578E-02  O.130O5E-O1  0.10401E+01 

0.80185E-02  0.15084E-01  0.10462E+01 

0.88649E-02  0.17343E-01  0.10524E+01 

0.97001E-02  0.19824E-01  0.10590E+01 

0.10648E-01  0.22591E-01  0.10665E+01 

0.12406E-01  0.25776E-01  0.10764E+01 

0.14699E-01  0.29533E-01  0.10885E+01 

0.17111E-01  0.34100E-01  0.11024E+01 

0.19935E-01  0.39973E-01  0.11198E+01 

0.20812E-01  0.41385E-01  0.11244E+01 

0.22640E-01  0.42975E-01  0.11312E+01 

0.22707E-01  0.44512E-01  0.11344E+01 

0.23777E-01  0.46259E-01  0.11401E+01 

0.24981E-01  0.48145E-01  0.11463E+01 

0.26155E-01  0.50175E-01  0.11527E+01 

0.27821E-01  0.52411E-01  0.11605E+01 

0.29190E-01  0.54810E-01  0.11680E+01 

0.30999E-01  0.57465E-01  0.11769E+01 

0.32988E-01  0.60403E-01  0.11868E+01 

0.35395E-01  0.63682E-01  0.11982E+01 

0.38177E-01  0.67306E-01  0.12110E+01 

0.41377E-01  0.71396E-01  0.12255E+01 

0.44855E-01  0.75962E-01  0.12416E+01 

0.48999E-01  0.81156E-01  0.12603E+01 

0.50986E-01  0.83984E-01  0.12699E+01 

0.54002E-01  0.87132E-01  0.12823E+01 

0.56569E-01  0.90420E-01  0.12940E+01 

0.59962E-01  0.94031E-01  0.13080E+01 

0.63492E-01  0.97932E-01  0.13228E+01 

0.67176E-01  0.10208  0.13385E+01 

0.70755E-01  0.10651  0.13545E+01 


0.76910E-01  0.11170 
0.82571E-01  0.11716 
0.88776E-01  0.12313 
0.96438E-01  0.12975 
0.10452  0.13708 
0.11504      0.14531 


0.13772E+01 
0.13995E+01 
0.14238E+01 
0.14524E+01 
0.14832E+01 
0.15207E+01 


• 


OAMA 


8.i 


XN 


u.72137E-u2    .uV 

O.13O2u£+01 

.ICO 

0.13016E+C1 

.  1 5  Cj 

0.13009E  +  '.  1 

•  2O0 

0.12998E+C1 

.250 

0.12982E+C1 

•  3i.i. 

o.l2958E+ol 

.35o 

0.12925E+-1 

•  400 

Q.12876E+U1 

•  45o 

0.128o6£+ol 

.500 

0.127O7E+L1 

PR 


63 

p 


7.5  0.84253E-O2 


.55o  0.12563E+01 
.600  0.12356E+01 
.650  0.12057E+C1 
•  7uO  0.11625E+01 
.750  O.K  999E  +  C  1 

.800  0.10092E+Q1 

.85o  0.87742 

.900  0.66578 

.95^  U.4V.678 

l.Oo^  0.12C7OE-03 


i.O  0.72137E-o2   l.OuO 


•  05o  0.1325OE+C1 
.100  U.13244E+01 
.150  0.13234E+01 
.200  0.13219E+01 
.250  0;l3197E+ol 

•  3uO  0.13166E  +  L  1 
.35u  0.13121E+01 

•  4oo  U.13059E+01 
.45v  0.12971E+01 
.500  0.12849E+01 

.550  o.12676E+^  1 
.600  D.12434E+01 
.650  0.12093E+! 1 
.700  0.11611E+01 
.750  0.1093UE+01 

•Sou  0.99666 
.85^-  0.66014 
.9ou  u. 66646 

•  95v.  0.39136 
l.OOo  0.2O27OE-C3 


7.5  u.84253E-o2  l.OOU 
7.0  0.99027E-O2  1.000 
6.8  0.10786E-U]  l.OOO 

6.6  0.11531E-O1  l.CLO 
6.4  0.12334E-ol  l.OuO 


j. 12714 


P. 15453 


0.1*120 

0. 15766 
0. 16727 
0.1753CJ 

0. 16  38; 


0.16461 
0.17627 
0.18180 
o. 18702 
0.19258 


0.15633E+01 


0.1612 OE+01 
0.16678E+01 

0.16981E+01 
0.17246E+;! 
0.17528E+01 


64 


KAMA 


XN 


px 


PR 


6.2  0.1327XE-01 

1.000 

6.0  0.14313E-01 

1.000 

5.8  0.15465E-01 

1.000 

5.6  0.16716E-01 

, 1.000 

5.4  0.18092E-01 

1.000 

5.2  0.19668E-01 

1.000 

5.0  0.21487E-01 

1.000 

4.8  0.23337E-01 

1.000 

4.6  0.25455E-01 

1.000 

4.4  0.28020E-01 

1.000 

4.2  0.30859E-01 

1.000 

4.0  0.33536E-01 

1.000 

3.8  0.36722E-01 

1.000 

3.6  0.40524E-01 

1.000 

3.4  0.44671E-01 

1.000 

3.2  0.49436E-01 

1.0^0 

3.0  0.54582E-01 

1.0J0 

2.9  0.57869E-01 

1.000 

2.8  0.60976E-01 

1.000 

2.7  0.63876E-01 

1.000 

2.6  0.67533E-01 

l.OuO 

2.5  0.71746E-01 

1.000 

2.4  0.75624E-01 

1.000 

2.3  0.79693E-01 

1.000 

2.2  G.84457E-01 

l.ooo 

2.1  0.88231E-01 

1.000 

2.0  0.95300E-01 

1.000 

1.9  0.97764E-01 

1.000 

1.8  0.10291 

1.000 

1.7  0.10916 

1.000 

1.6  0.12491 

1.000 

1.5  0.12578 

1.000 

1.4  0.14670 

1.000 

1.3  0.15288 

1.000 

1.2  0.14333 

1.000 

1.1  0.18319 

1.000 

1.0  0.20054 

1.000 

.9  0.16526 

1.000 

.8  0.18682E-01 

1.000 

.7  0.66941 

1.000 

.6  0.24838 

1.000 

.5  0.39536E+01 

1.000 

.4  0.53580E+01 

1.000 

.3-0.51193E+02 

1.000 

.2  0.64759E+02 

1.000 

.1-0.29826E+05 

1.000 

0.19354 
0.20420 
0.21581 
0.22822 
0.24166 

0.25680 
0.27407 
0.29138 
0.31090 
0.33424 

0.35982 
0.38361 
0.41146 
0.44429 
0.47971 

0.51998 

0.56305 
0.59039 
0.61616 
0.64010 

0.67015 
0.70470 
0.73642 
0.76959 
0.80832 

0.83890 
0.89608 
0.91595 
0.95724 


0.19851 
0.20474 
0.21139 
0.21832 
0.22558 

0.23333 
0.24183 
0.25019 
0.25926 
0.26939 

0.27983 
0.28972 
0.30061 
0.31236 
0.32430 

0.33678 

0.34942 
0.35659 
0.36304 
0.36939 

0.37603 
0.38359 
0.38994 
0.39696 
0.40344 

0.40895 
0.41794 
0.42204 
0.42836 


0.10073E+01  0.43383 

0.11339E+01  0.44424 
0.11409E+01  0.44822 
0.13090E+01  0.45571 
0.13588E+01  0.46040 
0.12823E+01  0.46442 

0.16016E+01  0.47137 
0.17411E+01  0.47698 
0.14584E+01  0.47732 
0.29524  0.48051 
0.55116E+01  0.49384 

0.21087E+01  0.48966 
0.32528E+02  0.50966 
0.44274E+02  0.49438 

-0.39828E+03  0.47521 
0.49466E+03  0.49740 

-0.15259E+06  0.19272 


0.17841E+01 
0.18179E+01 
0.18544E+01 
0.18931E+01 
0.19345E+01 

0.19303E+01 
0«203ieE+0l 
0.20831E+01 
0.21403E+01 
0.22073E+01 

0.22793E+01 
0.23466E+01 
0.24241E+01 
0.25133E+01 
0.26080E+01 

0.27135E+01 
0.28249E+01 
0.28940E+01 
0.29584E+01 
0.30190E+01 

0.30924E+01 
0.31766E+01 
0.32527E+01 
0.33331E+01 
0.34235E+01 

0.34957E+01 
0.36280E+01 
0.36760E+01 
0.37712E+01 
0.38823E+01 

0.41563E+01 
0.41783E+01 
0.45295E+01 
0.46385E+01 
0.44935E+01 

0.51459E+01 
0.54361E+01 
0.48715E+01 
0.25315E+01 
0.13011E+02 

0.61968E+01 
0.67074E+02 
0.90536E+02 

-0.79461E+03 
0.99171E+03 

-0.30517E+06 


65 


d-°°d 


66 
3.3.  Results  and  Discussion 

The  results  of  the  computed  entrance  length,  velocity  dis- 
tribution and  pressure  drop  in  the  entrance  region  by  this  method 
are  shown  with  sample  output  data  in  Table  3.5,  Table  3.10,  and 
Table  3.13  and  with  plotted  curves  in  Pig.  3.2  and  Fig.  3.5. 

Since  the  computational  procedure  is  more  complicated  than 
that  of  Schiller's  method,  we  separate  it  into  three  detailed 
FORTRAN  programs.   However,  the  whole  computation  may  be  put 
together  into  one  FORTRAN  program  as  described  in  section  3.2, 
part  (ill).   A  more  accurate  pressure  distribution  may  be  ob- 
tained If  we  divide  the  range  of  R  into  100  equl-distant  Inter- 
vals instead  of  20.   However,  when  we  compare  the  results  ob- 
tained in  each  case  we  find  that  at  y  =  20,  the  error  is  only 
1.8  x  10   ,  so  that  we  can  conclude  that  the  division  of  range 
of  R  in  20  intervals  provides  sufficient  precision  for  computing 
the  pressure  drops. 

In  the  computational  procedure  the  values  of  Iq(y)  and  I-i  (v) 
are  not  computed.   Although  the  equations  representing  In(v)  and 
I1(y)  are  available,  the  convergence  of  the  series  to  obtain 
these  values  is  too  slow  to  compute  these  with  the  IBM  1620 
computer.   Therefore,  the  values  of  Iq(y)  and  J^y)  are  read  in 
from  an  available  reference  [3]. 

The  results  of  the  present  work  are  compared  with  those  of 
Langhaar  for  some  values  of  y  as  shown  in  Table  3.14.   Except 
for  the  results  In  the  range  of  Y  <  1.2,  the  results  are  nearly 
the  same  for  1.2  <  y  <  100.   In  the  present  work,  for  the  region 
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of  0  <  y  <  1.2,  a  reasonable  result  could  not  be  obtained.   The 
reason  Is  unknown.   In  the  near  future,  by  using  an  IBM  360,  much 
faster  than  an  IBM  1620,  we  can  oompute  Iq(y)  and  I,(y)  for  the 
whole  range  concerned  and  we  may  be  able  to  obtain  more  accurate 
results  than  formerly. 
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CHAPTER  k.      SCHLICHTING'S  MATCHING  METHOD 

J+.l.  Description  of  the  Method  1 2,  10] 

In  this  chapter  we  shall  consider  the  case  of  flow  in  the 
inlet  length  of  a  straight  channel  with  flat  parallel  walls.   The 
flow  field  is  assumed  to  be  laminar,  steady  and  two-dimensional 
and  all  fluid  properties  are  considered  to  be  constant.  We  also 
assume  that  the  velocity  in  the  inlet  section  is  uniformly  dis- 
tributed over  its  width,  2a  (Fig.  ^.1).   The  governing  boundary 
layer  equation  and  continuity  equation  can  be  written  as 

.2 


u  |S  ♦  v  *-** 


f   dx  +  V  I    2  '  llJ 


^  +  ^  =   0  .  (2) 

ax  x  ay  v  ' 

A  system  of  co-ordinates  is  used  in  which  the  axis  of 
abscissas  coincides  with  that  of  the  channel  (Fig.  ^.l).   For  the 
expansion  in  the  upstream  direction  we  shall  measure  the  ordinate 
y  from  the  center-line  of  the  channel,  whereas  for  the  expansion 
in  the  downstream  direction  the  ordinate  y*  will  be  measured  from 
one  of  the  walls. 

In  this  method  the  flow  field  is  divided  into  two  sections, 
the  upstream  and  the  downstream.   In  the  upstream  section,  the 
boundary  layers  develop  in  a  way  similar  to  a  flow  over  a  flat 
plate  with  a  pressure  gradient  in  the  direotlon  of  flow.   The 
stream  function  and  center  line  velocity  are  expanded  in  a  power 
series,  as  in  the  case  of  the  flat  plate.   The  solution  is  found 
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after  the  coefficients  of  the  series  are  determined  by  substi- 
tuting the  assumed  series  Into  equation  (1).   For  the  downstream 
section,  the  velocity  Is  assumed  to  be  the  sum  of  the  fully 
developed  parabolic  velocity  distribution  and  a  deviation  veloc- 
ity which  approaches  zero  as  x  — *•  °o.   This  deviation  velocity  is 
obtained  by  a  series  expansion  in  the  upstream  direction.   Having 
obtained  both  solutions  in  the  form  of  a  series  expansion,  they 
are  joined  at  a  point  where  both  solutions  are  valid.   In  this 
way  an  approximate  description  of  the  flow  field  in  the  entire 
region  Is  obtained. 

For  the  upstream  solution,  the  center  line  velocity  Is 
assumed  to  be  a  series  in  6 . 

u0(x)   =  u(l  +  Kj€  +  K2£2  +  ...)  (3) 

where 


/  t/x 
/a2u 


and  K,  ,  K-,...  are  constants  to  be  determined. 

Substituting  the  expression  in  equation  (3)  Into  the  pres- 
sure gradient  at  the  center  line  in  the  flow  direction,  we  obtain 

-Jfi    -    -oTG?   -    4f  K*<4  +  2K2)€ 

+  3(^2  +  k3)&2  ♦...].       (*) 

It  Is  also  assumed  that  the  stream  function,  Iff  (x,y* ) ,  can  be 
expanded  into  a  series  in  £ 


where 


Then 
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4>(x,y")       =      5«[«0(»l)    ♦«2f1(1)   ♦    •••]  ^a> 

*-  yS-y  • 

u(x.y')      .    |jfj-    .     u[f-(H)   +€fj(»l)   ♦   ...]  .  (5) 

v(x.y')     =     -||    =     -|f[6f0  +  2€2f1+  362f2+   ...] 

+  %fMf6*efi  +  *2f2  *•••]•  <6> 

Substituting  the  expressions  in  equations  (4),  (5),  and  (6)  into 
equation  (1)  and  collecting  terras  in  powers  of  6 ,  we  obtain 
differential  equations  for  fQ,  f -,,...   The  equation  of  the  first 
approximation  (I.e.  the  zero  power  of  6)  is 

2fg«  +  f0fg   =   0  (7) 

with  the  boundary  conditions 

fQ  =  f  •   =  0     at     H  =  0  , 

f «  — *■   1  at     \  — »-  »  . 

This  Is  identical  with  the  differential  equation  of  flow  over  a 
flat  plate  given  and  solved  by  Blasius.   The  second  approximation 
(I.e.  the  first  power  of  6 )    is  determined  by  the  equation 

2f"*  *  f  f"  -  f'f  *  2f"f   -   -K  (8) 

1*01    0101       1  *  ' 


7^ 


which  la  a  linear  differential  equation  of  third  order.   The 
boundary  conditions  are 

f1     =  f  •   =  0     at     r\_  =  0  , 
f  ■   =  K  at    n.   —>   »  . 

These  differential  equations  can  be  solved  numerically. 
With  the  values  of  fQ,  f,,...  and  Kx,  K2,...  being  found,  the 
solution  of  the  upstream  section  is  complete. 

For  the  downstream  section,  we  assume 

u(x,y)  =  uf(y)  -  u*(x,y)  (9) 

where 

uf(y)   «  £  5(1  -  *£) 

Is  the  fully  developed  parabolic  velocity  distribution  and  u*  is 
the  deviation  velocity,  which  is  assumed  to  be  in  the  form 

J£  =  0le"Xl6  4>i(|)  .  (10) 

Substituting  the  expressions  In  equations  (9)  and  (10)  into 
equation  (1)  yields  the  following  differential  equation  forl^: 

Yiv*  Mi11  -f^K+Yi]  =   °  <"> 

with  K.    as  Its  eigenvalue.   The  boundary  conditions  are 
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^1=^J  =  0     at     y  =  0  . 

v|»i  =  v^i  =  °    at    y  =  ±a# 

V|l£  =  1  at     y  -  0  . 

The  last  boundary  condition  is  chosen  arbitrarily.   This  is 
possible  since  c.  Is  still  free.   Having  these  boundary  condi- 
tions, equation  (11)  can  be  solved  by  a  power  series  method. 
Assume 

SVa">   =  ^Va^ 

■  A0  ♦  Al(f>  +  A2(f)2  +  "•  •  (12) 


¥!<!>   =  _2,  n  V$ 


n-1 


n=l 


Al  ♦  2A2(a}  *  3A3(a)2  +  •••  •  <12a) 


V|,-{X,  =   s  n(n-l)An(|) 

n=<s 


n-2 


=   2A2  +  6A3(£)  +  12A4(|)2  ♦  ...  ,         (12b) 


Vj/J'Cf)   =   S  n(n-l)(n-2)An(f)n-3 

=  6A3  +  2^(f)  +  60A5(f)2  +  ...  ,        (12c) 
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V|>JT(|)   =   E  n(n-l)(n-2)(n-3)An(f)n-4 


n=4 
=   2^A4  ♦  120A5(^)  +  360A6(£)2  +  ...  .      (12d) 

By  substituting  the  expressions  in  equations  (12)  to  (12d)  into 
equation  (11),  the  recursion  formula  for  A  is  obtained. 

3*.-l  (n-2)(n+l)An  -  (n-l)(n+2)An+2 
An+^  =        2(n+l)(n+2)(n+3)(n+t)       ' 

The  values  of  A  and  \,  can  be  obtained  from  this  recursion 
formula  together  with  the  boundary  conditions.   Thus  the  devia- 
tion velocity  is  defined  if  c.  Is  known,  and  c1  can  be  determined 
in  the  process  of  matching  both  solutions. 

The  pressure  drop  in  the  upstream  section  is  obtained  from 
Eule r • s  equat Ion 

-  j         du. 
.  I  dp.  _     _0  (13) 

f   dx  -  u0  dx  •  KX3' 

The  boundary  layer  thickness  is  very  thin  in  this  section  and  the 
flow  is  mostly  irrotational  In  the  central  core.   Integrating 
equation  (13)  gives  the  pressure  drop  between  any  two  locations 
in  the  upstream  section.   On  the  other  hand,  the  pressure  drop  in 
the  downstream  section  involves  the  perturbation  velocity, 
u*(x,y).   Substituting  the  assumed  velocity  distribution  in  equa- 
tion (9),  into  the  boundary  layer  equation  (1)  and  noting  that 
the  pressure  varies  with  x  only,  we  obtain 
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Integrating  equation  (1^)  from  the  Joining  section,  g*  to  g  gives 
the  pressure  drop  in  the  downstream  section 

2        2 
-6(62  -  e*  )  ♦  yi^J'dXa  X   -e  x   ). 


-2  " 


For  the  complete  determination  of  the  velocity  distribu- 
tions, as  well  as  the  pressure  drops  in  the  entire  entrance 
region,  only  o,  and  6*,  the  location  of  the  matching  section, 
remain  to  be  determined.   This  is  done  when  the  upstream  and 
downstream  solutions  are  matched.   The  conditions  that  the  center 
line  velocities  of  both  Isolations  be  equal  and  that  the  pressure 
gradient  jjp  be  identical  at  the  matching  section  provide  two 
simultaneous  equations  for  c,  and  6*  as  follows: 

-\  e2 

l  +  Kj(S  ♦  K262  +  K3f3  +  K^t*  =  |  -  ^e  X        .         (15) 


jg   (1  +  Kjfr  +  K262  +  K3€3  ♦  K^X^  +  2K26  ♦  3K3f2  ■»■  ^fe3) 

=  3  +  ox*  V|)  J*(D  .      (16) 

Therefore  c,  and  6*  are  determined  from  equations  (15)  and  (16). 
The  solution  is  then  complete. 

The  detailed  derivation  of  the  equations  Is  given  in  Li's 
report  [  lj. 
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k.2.   Computational  Procedures 

In  this  method  the  velocity  distribution,  pressure  drops  and 
entrance  length  In  the  whole  entrance  region  can  be  represented  ■ 
numerically  as  soon  as  the  following  constants  and  functions  are 
found.   Those  are: 

(1)  The  functions  f_(1)  and  their  derivatives  in 
equations  (5)  and  (6). 

(2)  The  constants  K   in  equation  (3). 

(3)  The  eigenvalue  x,  in  equation  (11). 

CO   The  constant  c,  and  the  value  of  g  in  equations 
(15)  and  (16)  at  the  matching  point. 

The  values  of  f0(1)  and  its  derivatives,  K,,  \j,  c,  and  6 
have  been  obtained  in  the  present  work.   We  will  proceed  through 
the  computational  procedures  in  detail  as  follows: 

(1)   The  Function  f0C0  and  its  Derivatives 

Equation  (7), 

f0f"  «■  2fg'  =  0  (7) 

with  the  boundary  conditions  given  by  equation  (7a) 
f0   =   f  J   =   0     at      ^  =   0 


(7a) 


f  •  -*  1 


has  been  solved  by  Blaslus  [5]  using  a  series  solution,  by 
Toepf er  [  6  ]  using  the  Hunge-Kutta  integration  method,  and  by 
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Howarth  [  7 1  and  others.   Here,  the  solution  Is  obtained  by  the 
Bunge-Kutta  Integration  method. 

To  Integrate  equation  (7)  numerically,  another  boundary 
value  at  "\  =  0  is  required  In  order  to  start  the  integration. 
Thus  we  have  to  use  a  trial  and  error  method  assuming  many 
different  values  to  start  the  integration  until  we  meet  the  third 
boundary  condition  f '  — *■   1  as  1  — >   ».   This  difficulty  can  be 
avoided  if  we  set 

f0(>O   =  AF(-J)     where     5  =  Bl\ 

where  A  and  B  are  arbitrary  constants,  then 
fj(»l)   =  ABF'tf)  . 

T«W      =  AB2F»(5)  .  (17) 


fJ'Cl)   =  AB3F"'(-5)  . 


After  substituting  the  above  relations  into  equation  (7), 
the  transformed  differential  equation  becomes 

AB3F"'  +  i  A2B2FP"   =0  (18) 

with  the  corresponding  boundary  conditions 

F  =  F«   =  0     at     -£  =  0» 

(18a) 

F'  — '  AB  as     1   —  "  • 

Equation  (18)  can  be  made  identical  with  equation  (7)  if  we 
set  A  =  B.   Then  equation  (17)  becomes 
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2F»«  +  FF.i  m     0  (19) 

with  the  boundary  conditions 

F  =  F'   =  0     at     •$  =  0  , 

(19a) 

P»  — *■  \  as    -^  — *■   »  . 

A 

Now  equation  (19)  can  be  Integrated  by  giving  an  arbitrary 
third  boundary  value  at  -^  =  0.   Solving  F(5)  numerically,  we  can 
obtain  the  A  value  from  the  relation 


A  =  Jftzt  • 

Substituting  the  A  value  Into  equation  (17),  we  finally  get  the 
right  starting  value  from  the  relation 

f"(0)   =  A3F"(0)  .  (20) 

Thus,  in  order  to  obtain  the  value  of  fjl(0)  by  integrating 
equation  (19),  we  assume  F"(0)  =  1.  This  can  be  so  assumed  be- 
cause A  Is  an  arbitrary  value. 

The  Runge-Kutta  method!  ^  ]  is  applied  here  to  accomplish  the 
numerical  integration  for  equation  (19) »  We  will  describe  the 
key  point  of  this  method  as  follows 

H  =  ?•   =   f(5.F) 

with  Initial  values -J.  and  FQ. 

The  Increment  for  advancing  the  dependent  variable  is  given 
by 


81 


where 


AP  =  £  (Kj  ♦  2K2  +  2K-  ♦  K^) 


Kl  =  hf^o-  F0> 

K2  =   hf(-$Q  *  |  h.   PQ  +  |  Kl) 

K3  =   hf (^0  +  2  n'   F0  +  2  K2) 


K 


4 


=  hf(-§Q  ♦  h,   FQ  +  K3) 


The  values  at  ("£ii  F-i )  are  then  given  by 

-Si  =-?o  +  h'     Fl  =  F0  +  AF  • 

The  Increment  on  F  for  the  second  Intervals  Is  computed  by  the 
same  formulas,  with  (-§.,  FQ)  replaced  by  (■£,,  F,) .   Thus  all 
Intervals  are  computed  In  the  same  manner  by  using  for  the 
Initial  values  those  at  the  beginning  of  each  interval.   In  order 
to  facilitate  the  computation,  let  us  employ  a  double  subscript 
notation  in  which  the  first  subscript  denotes  the  interval  In 
which  we  are  working  and  the  second  the  entries  in  the  interval. 

-5n=*o'   Fii  =  Fo   «*  Fii  =  f<VFo>  =  f<*Lrpn>- 

Now  we  see  that 


and 


Ki  -  hf<V  V   -  wh 


Fo  *  I  Ki  =  Fo  *\  **h  =  po  +  T-  Fii  • 
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Thus  we  have 

5l2  ="^11  *  2  •   F12  =  Fll  +  ~Z    Fll'   F12  =  f*a2'  F12*' 
Continuing,  we  see  that 

K2  =  hf«0  +  f  h,   FQ  +|Kl)   =A"f  F'2 

K3  =  hf(50  +|h.   PQ  ♦  iK2)   =A-fP'3 

K^  =  hf(^0  *  h,   FQ  +  K3)   =  A-f  F'^ 

so  that 

AFi  =  ir(Fii  +  2Fi2*  2Fi3  +  FU> 

The  scheme  of  this  method  applied  to  solve  equation  (19)  is 
shown  partly  in  Table  4.1.   Using  this  scheme  it  is  easy  to  solve 
(19)  and  (7)  because  the  equations  are  the  same,  only  one  of  the 
boundary  conditions  (i.e.  Fw(0)  and  f"(0))  is  different.   A  flow 
chart  (Table  4.2)  is  then  constructed.  First  assume  F"(0)  =  1 
and  find  the  true  value  of  fg(0),  then  replace  F"(0)  by  f5(0)  and 
proceed  with  the  same  steps;  the  function  f0(1)  and  its  deriva- 
tives can  be  solved  with  no  difficulty.   The  program  for  finding 
f"(0)  and  for  solving  Blasius'  differential  equation  is  shown  in 
Table  4.4.   The  numerical  results  for  fg(0)  and  f^Cl)  with  Its 
derivatives  are  shown  in  Table  4.5. 
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Table  4.2.      Flow  chart  for  finding  fo(0)  and  solving 
Blasius     differentia!    equation. 


i  =  i 


I 


F»»l 


I 


Define 
r»i<  F«i 


Compute 
lit 
Fnr     by  using 


Compute 

F«».  Fn,',  F«», 
F«."  F.»,F«,' 
F«'i  .F.i^F'"* 

Fn^F-K^Fn* 


Compute 


Punch 

A,   fo'Vo; 


T 


I-I*i 


Fz  =  f.Co) 


>o 


Fm  *  Fm  *a  Fn 
Fiii  '  Fni  -i-^Fn 
r"i  *  Ftii  t^ET 


\, 

4 "(o)  - 

A3 

/' 

T»»"'  "■  n 

i 

• 

^1           ' 

J    Fn.-AF»i 

end}* — STopW 


>o 


<0 


w~* 


Pu/ic^ 

\i,  F^F*',  /F,*,J  Fi 


K=  kh 


_E 


«•  =o 


Table  ^.3.   Symbol  Table 
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Variable 

FORTRAN 
Symbol 

Variable 

FORTRAN 
Symbol 

-* 

E 

Pnk 

FC 

A-S 

DE 

Pnk 

F1C 

pm 

F 

Fnk 

F2C 

Fni 

Fl 

Fni 

F3C 

Fnl 

F2 

*Fn 

DF 

Fni 

F3 

A  F' 

DF1 

Fn2 

FA 

AFJ 

DF2 

Fn2 

F1A 

A 

A 

rn2 

F2A 

fg(0) 

F20 

Fn2 

F3A 

Pn3 

Ffl 

P» 
rn3 

FIB 

n3 

F2B 

F»3 

F3B 

Note:   When  1=2  in  the  Table  k.2, 
into  functions  f.  and  ■£  is 


all  functions  F  are  changed 
changed  into  i\. 
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TABLF  lt.it.       FORTRAN  PROGRAM  FOR  FINDING  f."(  0 )  AND  SOLVING 
BLASIUS  DIFFERENTIAL  EQUATION 

C   MATCHING  METHOD FINDING   F20   AND  SOLVING  BLASIUS  DIFF.  EOU. 

101  F0RMATI5E14.6) 

102  F0RMATI/3H  A=E14.6) 

103  F0RMATI/5H  F29=E12.6) 
1  =  1 

F2=1.0 

1  E  =  .0 
DF=.01 
F=.0 
F1  =  .0 

2  F?=-F*F2/2. 
FA=F+Fl»DF/2. 
F]A=F1+F2*0E/2. 
F2A=F2+F3*DE/2. 
F3A=-FA*F2A/2. 

3  FB=F+FlA*DE/2. 
FlB=Fl+F2A»DE/2. 
F2B=F2+F3A*DE/2. 
F3B=-FB*F2B/2. 

4  FC=F+F1B*DF 
F1C=F1+F2B*DF 
F2C=F2+F3B*DE 
F3C=-FC*F2C/2. 

5  DF=(Fl+2.*FlA+2.*FlB+FlC)»DE/6. 
DFl=(F2+2.*F2A+2.»F2B+F2C)»DE/6. 
DF2=(F3+2.*F3A+2.»F3B+F3C)»DE/6. 
PUNCH  101.  E.F.F1.F2.F3 

6  E=E+DC 
F=F+DF 
F1=F1+DF1 
F?=F?+DF2 

61  IF(DF1-.100000E-10)8.8,2 
B  K  =  0 
K  =  KT  +  1 

IFIK-1 )9.9.10 
9  A=SQRTF( 1./IF1-DF1 ) ) 
■  PUNCH  102,  A 
F20=A*A*A 
PUNCH  103.  F20 
1  =  1+1 
F2=F20 
GO  TO  1 
10  STOP 
FND 
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TABLF  4.5. 


NUMERICAL  RESULTS  CF  -f/u  )  AND  oOLUTICN 
BLA5IUS  DIFFERENTIAL  EQUATION 
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0.714331 
0.781206 
0.850563 
i  .922305 

- .996327 

o.  li  7252E  +  1  1 

o. 115076E+O1 

o.  1231ooE+ul 
0.131 305E+ol 

0. 139683E+ol 
0. 148222E+01 
0. 156912F+. ] 
0. 165740E+I 1 
0.174697E  +  . ,  1 


O.OOOOuO  0.332^63 

0.332061E-O1  0.332054 

0.664089E-01  0.331990 

C.9960O3E-O1  0.331815 

0.132766  0.331476 


0. 165888 
0.198941 
0.231894 
0.264714 
0".  297359 

0.329786 
o. 361944 
-.393783 
0.425244 
0.456269 

u. 480797 
0. 516765 
0.546110 
0.5  7  4766 
0.6  2 0/6 

0.629776 
0.656009 
0.681321 
0.7:  5664 
0.728994 

0.751 272 
0.77246  7 
0. 792556 
0.811522 
0.829356 

0.84605  7 
o  .  8  6 1  6  3  2 
876095 
889466 
9  177  5 


0.330917 
0.330085 

0.328928 
0.327395 
0.325439 

0. 323013 

0.32Cu77 
(  .316595 
...312  534 
0.307871 

0.302586 

0.296669 

0.29oll7 
J. 282950 
u. 275141 

0.266756 
0.257813 
o. 248355 
o.2 384 3 
B  .228095 

o.2 174 15 
v.2o6458 
0.195297 

O.184u09 
0.1 72671 

..'.161362 
o.15u157 
o. 139129 
(  .128348 
...11787  7 


F3 

0.000000 

-0.275656E-03 
-0. 110239E-02 
-C.247894E-02 
-0.440205E-02 

-0.686545E-02 
-0.985892E-02 
-0.133674E-01 
-0.173699E-01 
-0.218394t-01 

-0.267413E-01 
-0.320339E-01 
-0.376673E-01 
-0.435841E-01 
-0.497191E-01 

-0.560002E-01 
-O.62349ot-01 
-0.686817E-01 
-0.749110E-01 
-0.809472E-01 

-0.867003E-01 

-0.920821E-O1 

-0.970081E-01 

-0.101400 

-0.105187 

-0.  K'8308 
-0. 110715 
-0. 112372 

-0. 113257 
-0.113363 

-0. 112698 

-0.111283 
-0.109155 
-0. 1"6362 
-0.102964 
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0.35O00OE+O1 
0.360000E+01 
0.370000E+01 
0.38000OE+01 
0.39OO0OE+01 


0.183772r+01 
0.192955L+01 
0.202236E+01 
0.211606E+01 
0.221056E+01 


Fl 

0.913054 
0.923'43 
0.932O86 
0.941131 
0.948729 


F2 


F3 


0.400000E+01   0.230577E+01   0.955532 


0.410000E+01 
0.420000E+01 
0.430000E+01 
0.440000F+01 

0.450000E+01 
0.460000E+01 
0.4700Q0E+01 
0.480000E+01 
0.490000E+O1 

O.50COO0E+O1 
0.510COOE+01 
0.520000E+O1 
0.530000E+01 
0.540000E+01 

0.550000E+01 
0.560000E+01 
0.57.0000E  +  01 
0.58OOO0E+O1 
0.590000E+01 

0.600000E+01 
0.610000E+01 
0.620000E+01 
0.630000E+01 
0.640000E+01 


0.240164E+01  0.961594 

0.249807E+01  0.966970 

O.259501F+O1  0.971716 

0.269239E+01  0.975884 


0.279017E+01 
0.288828E+01 
0.298669E+01 
0.308536E+01 
0.318424E+01 

0.328331E+01 
0.338254E+01 
0.348191E+01 
0.358138E+01 
0.368096E+01 

0.378061E+01 
0.388033E+01 
0.398011E+01 
0.407992E+01 
0.417978EJ-01 

0.427967E+01 
0.437957E+01 
0.447950E+01 
0.457945E+01 
0.467940E+01 


0.979527 
0.982696 
0.985440 
0.987802 
0.989827 

0.991555 
0.993020 
0.994258 
0.995298 
0.996168 

0.996891 
0.997490 
0.997983 
0.998388 
0.998717 

0.998985 
0.999201 
0.999374 
0.999513 
0.999623 


0.650000E+01   0.477937E+01   0.999711 


0.660000E+01 
0.670000E+01 
0.680000E+O1 
0.690C00E+^1 

0.700000E+01 
0.710000E+U1 
0.720000E+01 
0.730000E+01 
0.740000E+01 


0.487934E+01  0.999780 

0.497932E+01  0.999833 

0.507931T+01  0.999875 

0.51793CL+01  0.999^08 


0.527929E+01 
0.537928E+01 
0.547928E+01 
0.557928E+01 
0.567927E+01 


0.999933 
0.999952 
0.999967 
0.999978 
0.999987 


0.107773  -0.990289E-01 

0.980870E-01  -0. 946319E-0 1 

0.888599E-01  -0. 898533E-01 

0.801263E-01  -0.847759E-01 

0.719119E-01  -0.794827E-01 

0.642342E-01  -0. 740549E-0 1 
0.571028E-01  -0.685701E-01 
0.505197E-01  -0.631009E-01 
0.444800E-01  -0. 5771 29E-01 
0.389725E-01  -0. 524646E-01 

0.339807E-01  -0.474059E-01 
0.294836E-01' -0.425784E-01 
0.254562E-01  -0. 380149E-01 
0.218710E-01  -0.337399E-01 
0.186981E-01  -0.297696E-01 

0.159066E-01  -0.26U32E-01 

0.134649E-01  -0. 227728E-01 

0.113416E-01  -0.197452E-01 

0.950575E-02  -0. 1 70219E-0 1 

0.792751E-02  -0. 145904E-0 1 

0.657845E-02  -0. 124353E-01 
0.543183E-02  -0. 1 053B7E-01 
0.446275E-02  -0. 888 1 12E-02 
0.364831E-02  -0. 744242E-02 
0.296766E-02  -0. 620208E-02 

0.240197E-02  -0. 5 1 3980E-02 

0.193442E-02  -0.423596E-02 

0.155012E-02  -0.347188E-02 

0.123597E-02  -0. 283003E-02 

0.980577E-03  -0. 229426E-02 

0.774077E-03  -0. 184980E-02 

0.608018E-03  -0. 148336E-02 

0.475201E-03  -0. 1 1 8309E-02 

0.369545E-03  -0. 9385 16E-03 

0.285947E-03  -0. 740503E-03 

0.220157E-03  -0. 58 1 1 37E-03 
0.168659E-03  -0.453632E-03 
0.128562E-03  -0. 35221 5E-03 
0.975096E-04  -0. 272016E-03 
0.735884E-04  -0. 208964E-03 
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E  F  Fl  F2  F3 

O.750O00E+C1  0.577927E+01  0.999993  0.552586E-04  -0. 1 59677E-03 

0.760000E+U1  0.587927E+01  0.999998  0.412875E-04  -0. 1 2 1 370E-03 

0.770000E+01  0.597927E+01  0.100000E+01  0.306949E-04  -0. 9 1 7666E-04 

0.780000E+01  0.6O7927E+01  0.100000E+01  0. 22706 1E-04  -0. 6901 82E-04 

0.790000F+01  0.617927E+01  0.100000E+0]  0.167127E-04  -0. 5 1 6362F-04 

0.800000E+01  0.627927E+01  0.100001E+01  0.122400E-04  -0. 384290E-04 

0.810000E+01  0.637927E+01  0.100001E+01  0.891954E-05  -0. 2R4501 E-04 

0.820000E+01  0.647927E+01  0.100001E+01  0.646744E-05  -0. 20952 1 E-04 

0.830COOE+01  0.657927E+01  0.100001E+01  0.466607E-05  -0. 1 53497E-04 

0.84000CE+01  0.667927E+01  0.100001E+01  0 . 334964E-05  -0. 1 1 1866E-04 

0.850000E+01  0.677927E+01  0.100001E+01  0.239262E-05  -0. 8 1 1 012E-05 

0.860000E+01  0.687927E+01  0.100001E+01  0.170051E-05  -0. 584913E-05 

0.870000E+01  0.697927E+01  0.100001E+01  0.120257E-05  -0.4 1 9654E-05 

0.880000E+01  0.707927E+01  0.1O0001E+01  0.846200E-06  -0. 299524E-05 

0.890000E+01  0.717927E+01  0.100001E+01  0.592465E-06  -0. 2 12673E-05 

0.900000E+01  0.727927E+01  0.100001E+0]  0.412745E-06  -0. 150224E-05 

0.910000E+01  0.737927E+01  0.100001E+01  0.286107E-06  -0. 1 05563E-05 

0.920000E+0)  0.747927E+0]  0.100001E+01  0.197335E-06  -0. 737961 E-06 

0.930000E+01  0.757927E''-01  0.1C0001E  +  01  0.135428E-06  -0. 5 1 3222E-06 

0.940000E+01  0.767927E+01  0.100001E+G1  0.9247B5E-07  -0. 355084E-06 

0.950000E+01  0.777927E+01  0.100001E+01  0.628350E-07  -0. 244405E-06 

0.960000E+01  0.787927E+01  0.100001E+01  0.424807E-07  -0. 167358E-06 

0.970000E+01  0.797927E*01  0.100001E+01  0.285765E-07  -0. 1 14010E-06 

0.98O00OE+O1  0.807927E+01  0.100001E+01  0.191274E-07  -0. 772678E-07 

0.990000E+01  0.817927E+0]  0.100001E+01  0.127389E-07  -0. 520974E-07 

0.100000F+02  0.827927F+01  0.100001E+01  0.844180E-08  -0. 349460E-07 

0.101000E+02  0.837927E+01  0.10C001E+01  0.556631E-08  -0. 233208E-07 

0.102000E+02  0.847927T+01  0.100001E+01  0.365198E-08  -0. 1 54831 E-07 

O.103000E+O2  0.857927L+01  0. lCO^OlE+01  0.238406E-08  -0. 1 02268E-07 

0.104000E  +  02  0.867927E  +  01-  0.100o01E+01  0.154859E-08  -0. 672030E-08 

0.105000E+02  0.877927E+01  0.100001E+01  0.100088E-08  -0.439350E-08 
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(2)   K-l  Value 

By  substituting  the  expressions  In  equations  C4-),  (5)t  and 
(6)  Into  equation  (1)  and  collecting  terms  In  powers  of  € ,  we 
obtain  differential  equations  (7)  and  (8) 

2fg'  +  f0f«  =  0  (7) 

2f»«  +  fQf»  -  f'f«  +  2fgf1  =  -^  .  (8) 

We  can  go  further  by  collecting  6J   terms  and  obtain  another 
differential  equation  such  as 

(fp2  +  2f'f£  -  fQf»  -  2fxfJ  -  3f2f5  "  2f|'   =  Kl  +  2K2     (21) 

and  so  on. 

The  boundary  conditions  for  these  equations  are  determined 
from  u=0  and  v=0  at  the  wall  (see  equations  (5)  and  (6)),  and 
become 

fQ(0)   =  f^O)   =  f2(0)   =  ...   =  fn(0)   =  0 

f'(0)   =  f«(0)   =  f£(0)   =   ...   =   fn(0)   =   0  . 

At  the  edge  of  the  boundary  layer,  H.~ »  and  u  =  u.(i)  so 
that  equating  equations  (3)  and  (5)  and  collecting  coefficients 
of  like  powers  of  €  yields 
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fJU)  =  l 
f{(«)   =  K1 


fj(-)   =  Kn  . 


For  large  \,    the  foregoing  relations  may  be  Integrated  to  yield 
the  approximate  relations 

f  0  ~  1  ♦  Ao 

f  1  ~  Kl^  +  Al 
f  0  ^  K„tt  «■  A, 


2 


K2^+  A2  (22) 


n      n  v    n 


Further,  from  the  continuity  equation  one  obtains 


=  /U  udy  =  [M»(x,y«)l 

0  V  = 


udy  =  !M»(x,y«) 

*y'=a 

At  the  center  line  of  the  duct  \   is  large  and  equations  (22)  are 
valid.   In  addition,  at  y=a  we  find  that  \  =  1/e  .   Using  these 
relations  in  the  foregoing  equation  together  with  equation  C+a), 
we  obtain 

K  K 

^U.y')   =  aa[€(i  +  AQ)  +  fe2(-^  +  A1)  +  63(-^  ♦  A2)  +  ...] 

1      =     [(1  +fA0)    ♦    (K^   ♦  6^)    +    (K2€2  +  63A2)    ♦    ...] 
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or  K,  =  -A.,  K0  =  -A,  ,  K-  =  -A„  ...  K  =  -A  ,.  Therefore,  when 
1     0   2     13     2      n     n-1 

equations  (7).  (8),  (21)  and  so  on  are  integrated  numerically, 
the  function  f  will  asymptotically  approach  the  expression 

fn  =  Kn^-Kn+1  '2» 

for  large  values  of  f\.      In  this  manner,  the  K  -values  needed  in 
each  equation  can  be  determined  from  the  solution  of  the  pre- 
ceding equation  of  the  set  (7),  (8),  (21)  and  so  on. 

Now  let  n=0,  then  from  equation  (3)  KQ  =  1,  and  using  the 
last  output  data  in  Table  k. 5,  one  obtains 

K,  =  10.5  -  0.877927  =  1.72073  . 

K,  =  1.72  was  computed  by  Schlichtlng  [2], 

(3)   The  Eigenvalue  j^ 

The  function  \j[.  (*)  is  assumed  in  series  form  as  shown  in 
equation  (12).   The  coefficient  A  can  be  expressed  by  a  recur- 
sion formula  as 

3».1[(n-2)(n+l)A  -  (n-l)(n+2)A   „) 
An+*  "  2(n+l)(n+2)(n+3)(n+4) ~  when  "i°   <24> 


_.   r(n-6)(n-3)An  ,    -  (n-2)  (n^A^   „-| 
An  =  2*1 1 n(n-I)(n-2)(n-3) —J"*™  n>4.(2Ua) 

From  the  boundary  condition  V|l,(0)  =  0,  we  find  from  equa- 
tions (12)  and  (2^)  that  AQ  and  all  constants  with  even  sub- 
scripts are  zero.   Fromty'(O)  =  1  and  equation  (12b),  A  is 
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found  to  be  one.   This  leaves  two  boundary  conditions,  vy.(+  1)  = 
VJIM+  1)  =0,  to  determine  the  eigenvalue  \,  and  A,.   Since  these 
boundary  conditions  are  in  the  form  of  an  infinite  series,  there 
is  no  simple  method  of  calculating  X.,  and  A.,.   A  trial  and  error 
method  is  then  used.   If  a  value  of  a.,  is  chosen,  then  values  of 
A  can  be  found  in  terms  of  A.,.   Fifty  terms  are  used  by  test  in 
calculating  the  summation  of  the  series  to  assume  a  small  error. 
Using  the  boundary  condition  ^, (+  1)  =  0,  A,  is  found.   These 
values  of  A.,  and  K-.    are  then  checked  with  the  last  remaining 
boundary  condition V^i (+  1)  =  0.   If  these  values  do  not  satisfy 
the  last  boundary  condition,  a  new  value  of  x.  Is  chosen  and  the 
process  is  repeated  until  the  correct  value  of  K-,    is  found. 

The  technique  of  performing  this  trial  and  error  method 
proceeds  using  the  following  steps: 

(I)  Given  a  wide  range  of  \.    (-130  to  130)  and  A, 
(-200  to  200),  for  each  \,  ,    we  assign  a  set  of  A., 
values  and  calculate  ^(j,  1)  and  Vj>{  (±  1)  by  using  equa- 
tion (2^a)  together  with  equations  (12)  and  (12a)  and 
plot  curves  as  shown  in  Fig.  k.Z,   Fig.  k.J,    and  Fig. 
'K'K  From  these  figures  we  can  predict  that  the  value 
of  A-  may  lie  somewhere  in  the  neighborhood  of  the 
origin.   We  then  take  a  smaller  range  of  A,  values  and 
proceed  to  the  next  step. 

(II)  Set  the  range  of  A^  from  -JO   to  +30,  and  assign 
several  values  of  Ju.  For  each  value  of  ju,  calculate 
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HJi(±  1)  andV|*J(  +  1)  as  before  and  trace  the  Inter- 
section points  ofVl>,(±  1)  and  *¥{(  +  1)  for  each  X,  as 
shown  in  Fig.  h.'j.      From  this  figure,  we  can  assure 
that  the  value  of  A.,  is  located  somewhere  between  10 
and  20,  and  that  the  value  of  A-  is  located  some- 
where between  -10  and  10. 

(ill)  From  the  conclusion  of  steps  (i)  and  (11),  we 
can  now  make  a  program  to  find  the  true  value  of  A,, 
which  lies  between  -10  and  10.   This  Is  done  by 
assigning  a  great  many  values  of  \,    and  applying  the 
Chords  method  [9]  (or  linear  interpolation  iterative 
formula)  (see  Fig.  4-.  6)  to  find  the  required  value 
of  \,    and  A.,  which  can  satisfy  the  boundary  condition 

Vf-^t  1)  =  0  andl^(  +  1)  =  0.  We  do  not  use  the 
Newton  iterative  formula  because  the  first  derivative 
of yx(±   1)  which  varies  with  k^   is  difficult  to  find. 

A  flow  chart  to  find  true  value  of  X-.  and  A-  is  shown  in 
Table  k.S-,    it  Is  followed  by  a  symbol  table  (Table  k.7)    and  a 
FORTRAN  program  (Table  k.Q) ,    and  some  output  data  are  shown  as 
in  Table  k.9.     From  these  output  data  we  can  see  that  K-.    =  18.81, 
A-,  =  -3.0103658  are  the  values  which  can  satisfy  the  boundary 
conditions  better  than  any  other  set  of  \.  and  A-  values. 
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4  Y.ia) 


«f.(A2) 


A3     - 


^(Ai)A2-».(Aa)Ai 
^i(A,)  -   <MAi) 


Fig.  4.6.     The   Chords    Method  (or  linear 
interpolation    iterative   formula) 


100 


Table  4.6.    Flowchart    for    finding  *,and 
A3   by   trial    and    error. 


7li=7\,*.ol 


I 


Y  *i*i7-f1 1 — (5WRri 


K»  0 


J*Jr  I 


T 


Aj  «A(j; 


J 


/v 


4B«/I62j 


AM' As 

AL  =  I 
XN  ■  3 


I 


m  ■=  M<) 


iZZE 


Aj  */)B 


Compute 


T 


•  Ai--AA 


-<\ 


XN»XN*2 


Compute 


AL*AH 
A*'  AH 


Oi»'ds  Method 


5.  *  5NW 


>o 


SNb)  « f . 


*■     K=Krl 


5/  =  ft 
4A*  /4a 


5**  £ 

40-  /I, 


Si'SMW 


S««)=  f, 
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Table  ^.7.   Symbol  Table  for  finding  J^ 
and  A,  by  trial  and  error 


Variable 

FOHTHAN  Symbol 

*■! 

DM 

An 

AN 

A3 

A3 

4*1 

SFN 

Yi 

SPIN 

An-2 

AM 

An-k 

AL 
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TABLE  4.8.   FORTRAN  PROGRAM  FOR  FINDING  A|  AND  A3 
BY  TRIAL  AND  ERROR 


C   MATCHING  METHOD — SOLVING  DAMDA  AND  A3  Br  IRIAL  AND  ERROR 
110  FCRMAT(F12.5,E16.8,2E14.6) 

118  F0RMATI6E12.5) 

119  F0RMAT1/28X.3E14.61 
DIMENSION  A(6)»SN(6) 
READ  118. ( A( J) .  J=l .2 ) 
DM=17.99 

1  DM=DM+.ul 

9  J=J+1 

IFIDM-18.86)  10.10,81 

10  A3=A(J) 
IFtJ-2] 101,102.102 

101  AA  =  A( 1 ) 
A3  =  AA 
GO  TO  11 

102  AR=A(2) 
A3  =  AB 

11  Aw=A3 
AL=1. 
XN  =  3. 
bFM=l ,+AM 
SF1M=1.+3.*AM 

12  XN=XN+2. 
AM=3.*DM# | (XN-6.)*r<N-3.)*AL-ixN-2.)*iAN-3.)*AM| 
AN=AN/(2.*XN*(XN-1. )*(XN-2. )*<XN-3. > ) 
SFN=SFM+AN 
3F1N=SF1M+XN*AN 
IFIXN-50.)  2u,30,30 

20  AL=AM 

AM  =  AN 

SFM=SFN 

SF1M=SF1N 

GO  TO  12 
30  <=K+i 

IF(K-,M31  ,32,34 
3  1  S  N (  i  ) -  '  F  H 

si=sn'(  :  i 

GO  TO  9 
32  SN<2)*SFN 

S2=SNl2 ) 

PUNCH  119,,  DM.SNC 1  I  »bN(2) 

Ip.(SN(  1  )*SN(  2  1-.0  )  33.1,1 
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33  A3M=(Sl*AtJ-S2*AA)/(Sl-S2) 
A3=A3K 

gc!  rc  11 

34  PUNCH  110,  DM.A3.SFN.SF1N 

IF( AQStSFNl-.oOOOOOl 1340.34 0,341 
34 v  1F<  Ar>MSFlN)-.OOCU,01>81»81  .1 
341  IF(SFN*S1-. 0)35.61.  36 

35  S2=SFN 
A3  =  A3  .. 

•  GO  TC  33 

36  IF(SFN*S2-.U)37,81,81 

37  51=5FN 
AA  =  A3 

GC  TC  33 
61  5TCP 
END 


input  data  cf  a<j> 

-  .10UCOE+02  0.1000UE+02 


TABLt    4.y.       SGMt    NUMERICAL    RESULT6    Oh    A|    AND    A3 
l>Y    TRIAL    AND    ERRSR 
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J.'-l 


6N<  1  ) 


SNI2J 


DM 


A3 


SHN 


SF1N 


•  181i>CE  +  02    -t,.2bl>3bblbt  +  Dl 


•181U00E+U2    -U.179290E+01 
O.K2330E-b7    -0*111614 


0.322494E+01 


•  16^utL  +  02  -u.2e /^2463t  +  vi 

•  182bCL  +  i'2  -o.2B /,524boE  +  i.l 

•  18200E  +  U2  -0.28  /i>2465E  +  Ul 


. 1B200UE+U2 

0.  12t  119t-Ub 

-0. Ib4294fc-0b 

-0.242940E-U  / 


-U.177211E  +  C11 
-U. 96138  ft-01 

-U.961409L-01 
-0.961402E-01 


0.32G240t+01 


•183G0E+O2  -0.28968282E+01 
•183oJt+02  -U.28968277E+0J 


.183C00E+O2 
-0.12763OE-U6 
-0.616171E-07 


-0.175145E  +  01 
-J.8O6B72E-01 
-'J.81.6683E-01 


0.318002E+01 


.1B4UGE+02    -u.2yla63olc  +  ,.i 
•  184^0t  +  vJ2    -u.  29iob^ociL+^i 


. 1S4U00E+02 

-0.  j41>4  /OL-ub 

0.9ea0i4L-i-o 


-0.173095E  +  01 
-G.651551E-01 

-U.b5ia44t-01 


0.315780E+01 


•  16bu^E  +  u2  -U.294^6523E+Ol 


0.4319'jbE-u  / 


-O.171o60c+01 
-u.49i>39bE-01 


O.313573E+01 


•1860CE+02  -0.296290UBE  +  (-l 
•187G0E+02  -0.29Bb37b2t+  01 


.  1860i^0E  +  o2 

-0.SS7723E-C C 

.  18  7COQE+02 

0.  /383<t6t-U7 


-0.  169u3oL  +  rjl 
-;;.3384S4E-01 
-0. 167032E+01 
-0.18069CE-01 


0.3U3B2E+O1 
0.309207E+01 


.1871CE+02  -u.298/b3fa^E  +  i.l 

•  lo/lut  +  i'2  -^.298 /63S6L+ul 

•18710C+02  -u.298 /b339L+0i 

•1871CE+02  -um citi  /6362E+01 


. 1871U0E+O2 
-0«32yl  /2L-^6 

0.  14CU4bt:-bO 

0. 1  /3G46E-^0 

-0.3479  WE-u  , 


-0.  X6bo32E  +  01 
-G.164b78t-Ql 
■0.1&46  72t-Oi 
•0  •  164o^  i'l-0  1 
-u.J.64bi>8t-01 


C.30B991E+01 


•  18720E  +  O2    -U.29B98991E+.J1 


.  lB72uOE  +  i-2 
0.  /b''62GE-b8 


-^.166632h  +  01 
-0.149G28E-01 


0.308774E+G1 


>1873ur:  +  02  -0  .  2^92  1  643  E+  u  1 
.18730E+U2  -u.29921632t+bl 


.  IB  /300t+02 
■U.2liV/  ^E-i/6 
-0.6b  7473E-0/ 


-0.166433L  +  '.;] 
-0.133198E-01 
-0.133206E-01 


(.'.308SS8E+01 


-io/'n-t+Oi  -^.2VV442VutTui 


.lb/4>.0E  +  u2 
o. a  1 J  uo4l-^ / 


-J.166233E+G1 
-u.  1 1  I  _>6uL-0l 


0.306342E+01 
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A3 


DM 
SFN 


jNtl  ) 

>F1N 


5«  ( 1 ) 


•  187!>k.GE  +  ^2  -O.166C34E+01       C.306125E+01 

.187bOE+G2    -O.29967006E+G1    -C.991 142E-G7  -i,.UH86E-01 

.1876G0E+G2  -U.165836E+01       0.307910E  +  01 

•  18760E+O2    -u.2V989716L  +  t>l       0.694117E-O8  -G .8 i6280E-02 

•187700E+U2  -G.16i636t+01       G.307694E+01 

.18770c+02    -u.3oul24boE+ul      0. 248Q1 7E-G6  -0.697654E-02 

.18770E+02    -u.3vul2461E  +  bl       u.  19i-i90E-w7  -u  .6W484E-02 


.18780E+O2    -0.30C35229E+O1 
.18780E+02    -U.3J035215E+G1 


.187800E+O2    -G.165437E+-01 

-C. 3570C0E-06    -o  .  5 38865E-02 

0.8786O3E-O7    -0. 538V22E-Q2 


.lb79uOE+u2  -G.1G5238E+01 
•18790E+02  -0.300i>e014E*Gl  -G.296336E-U6  -u«37<9993E-02 
•18790E+02    -u.Sj^ibUGit+ui      G.1314u4E-u7    - J. 3 ^9881E-02 


C.307478E+C1 


0.3G7263E+01 


•la8uuh+02  -C  .  3GuaG8  loE  +  ul 


• 1880oOL  +  k-2 
0.749293E-G I 


0.16SG40E+01 
0.22il65E-02 


0.30704UL+01 


.  188l0t  +  i~.2  -O.«3Ol036b8E+Gl 
.J881CE+J2  -O.301U3653E+0I 


.18820E+O2  -u.3~l26518E+~l 
•18820E+G2  -u. 3^  126bOBh  +  Gl 
•1882GE+U2  -U.3ul26512E+i>l 


.186 


.  1881tOE+k,2 
-0. 139789E-06 
-0.651258E-U7 

• 1882G0E+U2 

-0.254073E-O6 

0. 1 76383E-U6 

-0.  18361  /t-u6 


-G.164o41t+01 
-0.61'915CE-03 
-0.618435E-03 

-0. 164643E+01 
U.9725O0E-03 
Q.973849E-03 
.972249E-03 


^L,u^    -u.3U1265dlE+01    -0.836173E-O7      O.972249E-03 


.1883GE+C2    -U.30149389E  +  01 


.1884DE+02 

.186406  +  !  2 


-'   .3J1V2282E  +  01 
-0.31   1  /2293E+G1 


•1383G0E+U2 
-0.755189E-G7 

. 1884O0E+U2 

0.258549E-O6 

-0.776116E-07 


-O.16444bt  +  01 
U.256444E-02 

-0.164247E+C1 
G.4I5825E-02 
0.415768E-02 


•188S00E+U2  -0. 164u49E+0 1 

•18850E+02    -0.3U195229E+G1    -0.47G679E-O6  O.&74V26E-02 

.16fa'jUt  +  u2    -u.3ul9t>21uE  +  Gl       G.276i)bOE-u6  0.&7i^8bE-02 

.188506+^2    -0.3ol9S21?t+Gl    -0. 7 7 /869E-G7  0.574yl6fc-02 


.1886u6  +  (.2  -U.302181S3E  +  1/1 
•  18860E  +  U2  -i).30218175E  +  Gl 
.188606+0  2    -0.3021 8 168E+U1 


•  Io66u0h  +  v2  -u.  163ui)2L  +  &l 

J.486249fc-u6  0.734686E-02 

-o. 266742E-U6  O.734497L-02 

0.865253E-O7  U.734666E-02 


0.3C6833E+01 


0.306618E+01 


0. 3064U3E+01 


0.306188E+01 


0.305974E+01 


U.301>739t  +  01 
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(4)  C1   and  €  at  the  Matching  Section 

As  soon  as  \,  and  A,  are  found,  the  values  of  C,  and  6  at 
the  matching  section  can  be  obtained  by  using  equations  (15)  and 
(16).   Eliminating  0,6      from  both  equations  gives  a  seventh 
order  equation  in  €  as  follows  ri]. 

g(«)   =  Kj  ♦e[K^+2K2-6-l|)£'(l)3  ♦  e2  [3K1K2+3K3+2K14»J'(1)] 

+  63[1<-K1K3+2K2+Wj++2K2^^(1)]  +  £/4'[5K1K^5K2K3 

+  2K3(|»J'(1)]  +  e5[6K2Ki++3K3+2K^'(l)] 

+  ?K3K4f6  +  4kJ;€7  =   0  (25) 

The  K,  value  has  been  obtained  in  the  present  work  as  shown  in 
section  4.2-(2).   The  K,,  K-,  Kj,  values  are  taken  from  the  work 
of  Hoidt  and  Cess  [8], 

For  finding  6  in  equation  (25),  the  Newton-Haphson  iteration 
method  [  9  ]  is  applied  (see  Fig.  4.7)  because  g'  (6)  can  be  easily 
obtained. 

After  we  get  6 ,  C.  is  obtained  by  substituting  €  into  equa- 
tion (15).   A  flow  chart  (Table  4.10)  for  the  entire  computation 
is  followed  by  a  symbol  table  (Table  4.11)  and  a  FORTRAN  program 
(Table  4.12).   The  computation  is  also  carried  out  on  an  IBM  1620 
digital  computer.   The  output  data  are  shown  In  Table  4.13. 
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e.        e 


£v.r  ok- 


9  to 


Fig.  4.7.    Newton-Raphson's   Method. 
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Table  4.10.  Flow  chart  for  solving  C.and  I 


(STAlrt_» 


X.  .18.81 
Ptl.-.3o\oiiS 


Compute 


I  «  I'l 

1 *    6  ' 

iX 

Compute 

/"  r. 

1 

1     > 

f    Puncn     \           s    pur 

i 

(c,.  I-120T3 

Kj.  -2.236 

■ 

t  *  loo 
J  *° 

Table  4. 11.   Symbol  Table  for  Solving  C^  and  B 


109 


Variable 

FORTRAN  Symbol 

£k*l 

EX 

6k 

E 

g(6k)  or  g(€k+1) 

G 

S'(«k)  °r  8'(fk+1) 

Gl 

Kl 

XK1 

K2 

XK2 

K3 

XK3 

K^ 

xm 

*1 

DM 

A3 

A3 

kV-d) 

SF3N 

An 

AN 

e 

EL 

Cl 

Cl 

' 


TARLE  4.12.   FORTRAN  MROGRAC  FOR  SOLVING  CI  AND  £ 
AT  MATCHING  SECTION 
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u 


.'•'.ATCHlNb    METHCD- 
F0RMATI4E14.6) 


-SOLVING    EHSILON    MiNlb    Ci    AT    MATCnlNG    SECTION 


1L1    FCRKATI2EU.6) 

102    FCRMATI/3H    E=fcl4.6) 

li  3    ("0RMATI/4H    C1=E14.6) 

UM=18.B1 , 

A3=-3.-lu36b3 
11     Af/  =  A3 

AL=1  . 

X,'4=3. 

Sf3M=6.*AM 
•  12    XN-XN+2. 

AN=3.*0M* ( (XN-6. >*<XN-3. l*AL-<XN-2. )  * < XN-3. )*AM > 

AN  =  AN/  (2.KXN*-  (XN-1.  ;  •(>  ,-2.  )«<XN-3  .  )  I 

bf-3(M  =  Sr  3M+XN*  (XN-1 • )  *<XN-2.  )  *A,N 

1  F  (  XN-&G.  }'2ii  »22  .22 
21;    AL  =  AM 

AM  =  AN 

SF3M»5F3M 

GO    TO    12 

22  XK1 "1.72073 
XK?=-?.?38 
X<?=]9.81 
X<4=-2L8.896 

23  ol=XKl*XKl+2.*X<2-6.-SF3N 
o2=3.*XKl*XK2+3.*Xk 3+2.*X<l*SK3N 
B3=<f.*XKl*XK3+2.*XK2*XK2+4.*XK4+2.*Xf.2*SF3N 
B4=!>.»XKl*XH.4  +  b.*XK2*XK.3+2.*XK.3*SF3N 

ob  =  b.~X<2*XK.4  +  3.*XK3*XK3  +  2.  »XM»*SF  3N 

B6»=7.*XK3*XK4 

u7=4.*XK4*XK4 

24  E=100. 
I=! 

3 0     I=i+1 

G»XK  l+E*8H-E*£*B2+E»fc*E*B3  +  E*E*E*E*b4+E*E*E*E*E*B5 

G=G+E*t*E*E*E*£*B6+E*E*E*£*F*E*E*b7 

G 1=01+2 •*B2*E+3.*b3*t*E+4.*B4*fc*f_*t+5.*D5*E*E*E*E 

bl=ul+6.sc)6*t*t«E*h#L  +  7.*b7*t*L*t:'"';i*t*E 
tX=t-G/Ul 

IF(ABS(G)-.lGL00E-04)    50 » 50 .35 

3d    IK  1-1  )     4-.4u,41 
40    PUNCH    Kui    XK1»XK2»XK3.X<4 
GO    TO    42 


Ill 


M  PUNCH  10] ,  E,G 
4  2  h-iX 

60  TO  30 
50  LL=2. 7182818 

P=DM*E»E 

0»l.+XKl*E+XK2*E*e+XK3*E*E«E+XK4*E*t*E*fc 
Cl= ( 1.5-Q)«EL**P 
PUNCH  102,  £ 
PUNCH  103.  CI 
80  STwP 


kNU 
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TABLET  4.13. 


SOME  OUTPUT  DATA  FOR  CI  AND  £ 
AT  MATCHING  SECTION 


XK.1 


XK2 


XK3 


•172u73L+^l    -i..2238k<ut  +  ul       u.lvaiuot  +  o^    -u.2ot>B96t  +  03 


.539884E+02 
.462792E+02 

•396713E+U2 
. 340374E+02 
•291526F+-2 
.2<49913E+*2 
.214245E+v2 
.183673t+v2 
. l£>7468E+w2 
.1351  i  7E+U2 
.1  15755E+V2 

•  992527F  +  </l 
.851082E+01 
.72984&E  +  --1 
.62b929E+ul 
.53686uE+ui 
.46u517E+ul 
•39b083E+-l 

•  339O0OE+0  1 
.29o933E+C-l 
.249738-E+uJ 

•  2I443i>fc+ul 
.1841S4E+01 
.  lb8267t  +  ol 
.  136f'68E  +  ol 
.117063E+O1 
.1O08O3E+O1 
.869(81 

•  7.5^592 
.649912 
•b6492o 
.494040 
.43632  J 
.391598 

•  36i.  582 
.344f  82 
.339413 
.339072 


0.232642E+18 
U.79079OE+17 
0.2688O3E+17 
0.913705E+16 
:..3105H4E+16 
0.  105572E+16 
0. 3bbdboE+lb 
O. 121982E+ 1 : 
0.414636E+14 
0.  J  40941E+14 
G.479C  82F  +  13 
0. 162843E+13 
0.S53543E+12 
.  1  a  8  1  b  I E  +  1  2 
0.639572E+11 
0.217398E+11 
G.738958E+10 
0.251177E+10 
0.853755E+09 
0.290187E+09 
0.986299E+U8 
U.33521uE+08 
0. I13918E+08 
0.38  /089t  +  <-7 
0.  13J504E+O7 
0.446601  E  +  L'6 
0. 1M5B3E  +  06 
0.51398i>E+O!> 
0.1/398uE+>;5 
*-.b870ti2E  +  i.4 
0.  196961E  +  L/4 
0.61>33b8E+03 
0.21168bfc+u3 
0.650096E+02 
v..  1  /434uE  +  i,2 
0.320316E+I 1 
0.204601 
<  .  1028C0E-I  2 


P=       (.'.339u/>, 


cl  = 


.  1881-16E  +  1.2 
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^,3.  Results  and  Discussion 

In  1912  Topfer  [6]  calculated  the  fgCl)  function  by  using 
the  Hunge-Kutta  method;  the  results  were  then  improved  with  a 
high  degree  of  accuracy  by  Howarth  [  7 3  in  1938.   In  the  original 
paper  by  Howarth,  the  method  used  in  the  computation  was  not 
mentioned.   By  comparing  the  present  results  with  those  of 
Howarth  we  find  that  there  is  excellent  agreement. 

The  computation  of  \,  and  A,  by  trial  and  error  is  extremely 
difficult;  one  has  to  make  a  great  many  tests  before  getting 
favorable  results.   Schlichting  [  lo]  was  the  first  to  find  the 
eigenvalue  \-,  =  18.75  in  1934-;  other  authors  also  use  the  same 
value,  as  can  be  seen  from  the  literature  [  3.  8],   But  from  the 
output  data  in  Table  4.9.  we  can  conclude  that  k,    =   18.81  is  a 
more  accurate  one  that  can  satisfy  the  boundary  conditions  more 
accurately. 

The  values  of  C,  and  €  at  the  matching  section  are  obtained 
from  the  \,  value,  so  that  the  values  of  C,  and  6  computed  in 
the  present  work  are  different  from  those  obtained  by  other 
authors  as  shown  in  Table  4.14.   Some  of  the  results  computed  In 
the  present  work  are  compared  with  those  obtained  by  other 
authors  as  shown  in  Table  1.14. 
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The  study  of  the  developing  velocity  profiles  and  pressure 
drops  In  the  hydrodynamlc  entrance  region  of  tubes  and  ducts  has 
attracted  the  Interests  of  investigators  for  many  years.   Al- 
though the  exact  analytical  solution  Is  not  available,  many 
approximate  solutions  have  been  presented.   The  approximation 
methods  of  solution  can  be  classified  into  four  different  cate- 
gories:  the  momentum  integral  method,  the  linearizing  method, 
the  matching  method,  and  the  finite  difference  method. 

The  purpose  of  this  report  is  to  present  a  detailed  compu- 
tational study  of  the  momentum  integral  method,  the  linearizing 
method,  and  the  matching  method  by  using  a  high  speed  digital 
computer.   The  computational  procedures,  flow  charts,  FORTRAN 
programs  and  results  in  each  method  are  presented  in  detail  and 
some  of  the  results  are  compared  with  those  obtained  by  other 
authors. 

The  main  numerical  methods  used  in  this  work  are  numerical 
integration,  series  solution  of  differential  equations,  eigen- 
value solutions  and  the  evaluation  of  modified  Bessel  functions. 
In  Chapter  2  Schiller's  solution  by  the  momentum  integral  method 
for  a  flow  in  a  circular  tube  Is  described.   It  is  then  followed 
by  computational  procedures  to  compute  the  velocity  profiles, 
pressure  drops  and  entrance  length  in  the  whole  entrance  region. 
All  computations  are  done  using  only  one  FORTRAN  program  and  are 
carried  out  on  an  IBM  1620  digital  computer. 

Chapter  3  contains  a  description  of  Langhaar's  linearization 
method  for  the  flow  in  a  circular  tube.   The  nonlinear  partial 


differential  equation  of  motion  for  the  axial  direction  is 
linearized  by  the  assumption  that  the  convectlve  terms  are  a 
function  of  axial  direction  only.   The  solution  Is  presented  in 
closed  form  In  terms  of  modified  Bessel  functions.   The  evalua- 
tion of  modified  Bessel  funotions  and  the  trapezoidal  rule  for 
numerical  integration  are  used  In  this  chapter  to  obtain  the 
numerical  results.   In  order  to  make  the  computation  clear,  three 
separate  FORTRAN  programs  are  used  to  compute  entrance  length, 
velocity  profiles  and  pressure  drops.   The  results  are  compared 
with  those  obtained  by  Langhaar. 

Chapter  ^  deals  with  Schlichting's  matching  method.   In  this 
method  the  flow  field  is  divided  into  two  sections:  the  upstream 
and  downstream  sections.   In  the  upstream  section,  a  boundary 
layer  formulation  of  equations  is  used  and  a  solution  developed 
in  a  series  stream  function  with  Blaslus*  functions  as  coeffi- 
cients.  When  this  solution  becomes  unwieldy,  an  exponential 
velocity  deviation  from  the  fully  developed  flow  Is  assumed  and 
Joined  to  the  boundary  layer  solution  to  complete  the  description 
of  the  flow.   The  Runge-Kutta  integration  method,  the  Newton- 
Raphson  iteration  method,  and  the  Chords  method  (or  linear  inter- 
polation Iterative  formula)  are  applied  in  this  chapter  to  obtain 
the  numerical  results  for  the  entire  entrance  region.   The  re- 
sults are  then  compared  with  those  obtained  by  other  authors. 


